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Abstract 


These  notes  are  intended  to  serve  as  a brief  instruction  manual  for  engineers  engaged 
in  the  spectral  analysis  of  experimentally-derived  random  data.  This  manual  is  intended  to 
accompany  a software  routine  (Spectrum  Calculation  Routine,  V 1.0,  David  E.  Hess) 
designed  to  compute  amplitude  density,  autospectral  (power)  density  and  cross-spectral 
density  functions  for  discrete,  stationary  random  time  series.  This  program  will  run  on  any 
PC-compatible  desktop  computer  and  accepts  input  data  files  obtained  from  experiments  in 
which  analog  data  has  been  digitized  by  an  analog-to-digital  converter. 

A collection  of  techniques  are  described  which  are  necessary  for  the  proper  software 
implementation  of  spectral  analysis  procedures:  data  transformation,  mean  value  and  linear 
trend  removal,  digital  filtering,  time  history  tapering  and  data  overlapping.  The  various 
spectral  density  functions  are  then  defined  and  explained  along  with  a means  for  the  cal- 
culation of  the  error  associated  with  a particular  estimate.  Several  detailed  examples  are 
included  which  serve  to  illustrate  the  methods  described.  The  manual  also  contains  a brief 
description  of  the  considerations  necessary  when  sampling  the  data  to  be  analyzed.  Ref- 
erences to  more  complete  sources  of  information  are  provided  as  required. 


Notice  to  Users 

This  program  was  developed  for  use  with  the  ongoing  Compliant  Surface  Program  in 
the  Fluid  Flow  Group;  this  group  is  a part  of  the  Process  Measurements  Division  at  the 
National  Institute  of  Standards  and  Technology.  With  the  evolution  of  the  SPECTRUM 
routine  to  its  present  form,  the  author  felt  that  the  program  may  perhaps  be  of  some  utility 
to  others.  Therefore,  SPECTRUM  is  being  released  as  freeware.  This  program,  the  source 
code  and  accompanying  utilities  may  be  freely  copied  and  distributed,  but  may  not  be  sold. 
There  is  no  warranty  of  SPECTRUM’S  suitability  for  any  purpose,  nor  any  acceptance  of 
liability,  express  or  implied. 

Originally,  the  SPECTRUM  routine  was  designed  to  run  on  any  IBM-PC  or  compatible 
computer  using  the  Intel  80x86  series  of  microprocessors.  However,  the  code  was  compiled 
using  the  Microsoft  Fortran  Optimizing  Compiler,  V 5.00,  and  this  compiler  allows  the 
selection  of  a variety  of  options  at  compile-time  which  enhance  the  speed  of  operation  of 
the  program  but  limit  its  suitability  for  use  with  certain  computers.  The  version  contained 
on  the  diskette  requires  a computer  equipped  with  an  80286  or  better  microprocessor 
and  an  accompanying  80287  or  better  coprocessor.  A version  of  the  routine  which  will 
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run  on  computers  with  an  8086  or  8088  processor,  or  a version  which  does  not  require  a 
numeric  coprocessor  can  be  obtained  by  contacting  the  author.  The  program  uses  about 
450K  bytes  of  available  RAM,  and  requires  that  a hard  disk  or  RAM  disk  be  available  for 
storage  of  temporary  files.  The  amount  of  space  to  be  reserved  depends  upon  the  size  of 
the  data  file  to  be  processed,  but  can  easily  be  as  much  as  several  megabytes. 

Microsoft  has  announced  the  release  of  V 5.1  of  their  Fortran  compiler.  To  this  version 
the  capability  of  porting  programs  to  the  Windows  environment  has  been  added.  This  is 
significant  only  in  that  it  allows  the  use  of  the  DOS  extender  which  is  a part  of  the  Windows 
environment.  Fortran  programs  written  for  Windows  are  not  limited  by  the  DOS  memory 
restriction  of  640K,  but  instead  can  access  substantially  larger  amounts  of  available  RAM 
at  run-time.  This  can  greatly  increase  the  utility  and  speed  of  the  SPECTRUM  routine,  and 
a version  incorporating  these  improvements  may  be  available  at  a later  date. 

This  program  has  been  subjected  to  extensive  testing  and  appears  to  be  free  of  any 
bugs;  however,  it  is  possible  that  not  all  situations  have  been  anticipated.  The  author 
encourages  any  modifications  to  or  customization  of  the  source  code  as  required  by  the  user, 
and  welcomes  any  suggested  improvements  or  bug-fixes.  Please  direct  all  correspondence 
to  the  address  given  below. 


National  Institute  of  Standards 
and  Technology 
Building  230  Room  105 
Gaithersburg,  MD  20899 
Attn  : David  E.  Hess 

Phone  : (301)  975-5937 
Fax  : (301)  258-9201 
Email  : HESS@ENH.NIST.GOV 
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1 Overview  and  Parameters 


This  program  is  designed  to  compute  amplitude  density,  autospectral  (power)  density 
and  cross-spectral  density  functions  for  discrete,  stationary  random  time  series  in  the  form 

+ for  n = 0, 1,2, (1.1) 


where  N is  the  number  of  data  points  per  record.  Typically,  M such  records  are  stored 
in  a data  file,  and  each  record  is  a contiguous  segment  of  the  stationary  process.  Most  often, 
this  data  is  in  the  form  of  2-byte  integers  after  having  been  digitized  by  an  analog  to  digital 
converter.  This  routine  reads  this  data  file  and  then  transforms  the  integer  data  back  into 
voltages  or  optionally  into  some  other  desired  quantity.  This  transformation  into  a 
dimensional  quantity  is  by  means  of  a polynomial  of  up  to  fourth  degree.  The  coefficients 
necessary  to  complete  the  transformation  are  placed  in  coefficient  data  files  which  are  read 
by  the  routine.  Alternatively,  the  routine  can  read  4-byte  floating-point  data  which  has 
already  been  transformed  by  some  other  program.  The  input  data  may  consist  of  one  or  two 
separate  time  series  which  are  to  be  analyzed. 

After  transformation,  various  cosmetic  operations  may  be  performed  on  the  data  prior 
to  the  spectrum  computation.  The  mean  value,  which  is  calculated  during  the  transformation 
process,  may  be  subtracted  from  the  data;  or  instead,  a linear  trend  in  the  data  may  be  removed 
using  a least  squares  fit  technique.  Recursive  digital  filters  including  lowpass,  highpass, 
bandpass  and  bandstop  have  been  incorporated  and  may  be  applied  prior  to  the  calculation 
of  the  spectra.  Finally,  to  reduce  the  leakage  of  power  from  other  frequencies  due  to  the 
fact  that  the  input  time  series  are  of  finite  length  and  are  then  truncated,  a tapering  operation 
has  been  included.  A choice  of  four  tapering  functions  is  supplied,  and  the  appropriate  scale 
factors  are  computed  automatically.  The  use  of  tapering  may  introduce  some  additional 
variability  in  the  resulting  spectral  estimate;  to  counteract  this  tendency,  a means  for 
overlapping  the  input  data  records  is  provided.  This  option  employs  50%  overlap  to  reduce 
overall  error  at  the  expense  of  computation  time  (which  is  doubled);  however,  the  use  of 
overlapping  requires  that  the  data  be  contiguous  in  time  from  record  to  record. 

After  computation  of  the  selected  spectral  density  function,  the  data  are  stored  to  an 
ASCn  output  file  which  may  then  be  imported  into  a plotting  program.  The  routine  is 
designed  for  multiple  file  processing;  many  input  data  files  may  be  processed  sequentially 
if  certain  naming  conventions  are  followed.  The  maximum  number  of  files  is  currently  set 
to  IFMAX  = 100  in  order  to  allocate  storage  for  data.  Entering  data  file  names  on  the 
command  line  upon  invocation  of  the  program  is  also  permitted  and  the  maximum  is  set  to 
JFMAX  = 5.  Each  record  of  the  input  data  file  may  contain  no  more  than  NMAX  = 16384 
points;  note  that  in  the  case  of  two  channels  of  data  per  record,  each  channel  may  contain 
no  more  than  NMAX  / 2 points.  As  described  below,  the  order  of  the  recursive  digital  filters 
is  set  by  the  number  of  cascaded  stages  desired;  the  maximum  number  of  stages  is  currently 
placed  at  NSMAX  = 20.  For  convenience,  the  base  10  logarithms  of  certain  data  in  the 


output  are  computed;  to  avoid  taking  the  logarithm  of  zero,  any  number  in  the  output  data 
less  than  e - 1.0  x 10"^  is  set  equal  to  e . Any  of  these  maximum  values  may  be  altered 
by  making  a change  in  the  appropriate  PARAMETER  statement  near  the  beginning  of  the 
code  and  then  recompiling. 

The  speed  of  the  routine  is  limited  by  the  large  number  of  I/O  operations  that  are 
required;  reduce  the  number  of  points  per  record  and/or  the  number  of  records  per  file  for 
increased  speed.  Alternatively,  specify  in  the  configuration  file  that  the  routine  use  a RAM 
disk  for  storage  of  temporary  files.  The  size  of  the  executable  file  is  largely  governed  by 
the  setting  of  NMAX;  since  this  parameter  must  be  a power  of  two,  the  largest  possible 
value  is  16384  for  a DOS-based  computer  which  results  in  an  executable  file  size  of  430K 
bytes.  The  CHKDSK  command  must  report  that  at  least  450K  bytes  of  memory  are  present 
for  the  routine  to  execute  for  this  value  of  NMAX.  When  compiling  the  source  code  using 
Microsoft  Fortran  V 5.00,  at  least  603K  of  RAM  must  be  available  for  full  optimization  of 
the  code. 

When  the  data  to  be  processed  by  the  spectrum  routine  are  samples  of  a continuous 
function  of  time  at  a single  spatial  location,  then  the  calculated  spectral  density  functions 
are  a function  of  frequency.  The  user  should  keep  in  mind  that  if  the  input  data  are  samples 
of  a continuous  function  of  one  spatial  direction  at  an  instant  of  time  (such  as  intensity  data 
for  horizontal  or  vertical  rows  of  pixels  on  a photograph  firom  a CCD  camera),  then  this 
routine  can  be  used  to  determine  one-dimensional  estimates  of  wavenumber  spectra. 


2 Input  Data  Files 

For  any  routine  in  which  the  processing  of  large  amounts  of  data  is  required,  a rather 
rigid  structure  is  necessarily  imposed  on  the  manner  in  which  the  data  is  transferred  into  the 
program  for  manipulation.  The  spectrum  routine  is  certainly  no  exception  to  this  rule. 
Furthermore,  since  the  spectrum  routine  has  the  ability  to  process  data  files  in  batches, 
additional  rules  are  required,  such  as  naming  conventions,  to  allow  this  feature  to  be 
implemented.  Although  the  proliferation  of  rules  can  be  a bit  wearying,  later  sections  discuss 
some  situations  in  which  the  rules  can  be  bent  or  discarded. 

2.1  Naming  Convention 

A convention  for  naming  the  input  data  files  was  chosen  to  facilitate  multiple  file 
processing.  The  selected  method  resulted  fi-om  the  ease  with  which  the  file  names  could  be 
generated  by  a few  lines  of  code.  Each  file  name  consists  of  eight  characters:  the  first 
character  must  be  alphabetic  (A-Z),  the  next  three  characters  are  digits  (0-9)  and  the 
remaining  four  characters  are  a period  followed  by  the  letters  DAT.  Note  that  the  routine 
is  not  case-sensitive,  lower-case  as  well  as  upper-case  characters  may  be  entered;  however, 
all  lower-case  inputs  are  converted  to  upper-case  internally.  The  routine  automatically 
appends  the  suffix  [.DAT]  and  therefore  it  should  never  be  entered  upon  input.  The  data 
file  must  reside  in  the  same  directory  as  the  executable  file  unless  specified  otherwise  in  the 
configuration  file  (discussed  later).  Examples  of  acceptable  data  file  names  are:  A001.DAT, 
z087.dat  and  M418.DAT.  The  latter  filename  would  require  that  IFMAX  s 418.  Every 
rule  has  an  exception;  data  file  names  entered  on  the  command  line  when  calling  the  routine 
need  not  be  so  restrictive  as  outlined  above.  They  must  still  consist  of  eight  characters  and 
the  first  character  must  still  be  alphabetic  (A-Z);  however,  the  next  three  characters  can  be 
any  of  the  alphanumeric  characters.  The  blank  character  is  used  to  delimit  the  names  on 
the  command  line  and  is  therefore  not  a valid  character  for  a file  name.  The  remaining  four 
characters  must  be  a period  followed  by  the  letters  DAT.  Examples  of  acceptable  filenames 
for  files  entered  on  the  command  line  are:  TEST.DAT,  x_in.dat  and  S1$$.DAT.  When 
actually  entering  the  names  on  the  command  line,  only  the  first  four  characters  are  entered 
and  the  program  appends  [.DAT]. 

2.2  File  Structure 

All  input  data  files  must  be  FORTRAN  data  files  opened  as  sequential,  unformatted 
files.  The  minimal  amount  of  internal  formatting  allows  the  files  to  be  more  compact  which 
is  a consideration  when  dealing  with  large  quantities  of  data. 

2.2.1  File  Header 

The  first  record  of  the  data  file  consists  of  a header  containing  four  quantities  which 
characterize  the  data  which  follows.  The  first  quantity  is  a 2-byte  integer  which  gives  the 
number  of  channels  of  data  per  record  (1  or  2).  The  next  number  is  the  record  size  in  bytes 
and  may  be  either  a 2-byte  or  a 4-byte  integer.  The  record  size  is  twice  the  number  of  points 
per  record  if  the  data  consists  of  2-byte  integer  values,  or  the  record  size  is  four  times  the 


number  of  points  per  record  if  the  data  consists  of  4-byte  real  data.  The  third  quantity  is  a 
2-byte  integer  which  gives  the  number  of  records  of  data  which  follow,  and  the  last  quantity 
is  the  time  interval  in  microseconds  between  data  points  of  the  same  channel  expressed  as 
a 2-byte  integer.  Recall  that  the  largest  positive  signed  number  that  may  be  stored  as  a 
2-byte  integer  is  32,767  and,  for  a 4-byte  integer,  the  value  is  2,147,483,647. 

For  example,  suppose  the  data  file  consists  of  100  records  of  one  channel  of  data  saved 
as  2-byte  integers  and  sampled  at  10000  samples  per  second  with  2048  data  points  per  record. 
The  file  header  would  contain  the  four  values: 

1 4096  100  100 

Now  suppose  that  two  channels  of  data  were  sampled  with  a time  interval  of  0.016383 
seconds  between  each  successive  data  point  resulting  in  a time  interval  of  0.032766  seconds 
between  each  successive  data  point  of  the  same  channel.  Furthermore,  40  records  of  data 
were  sampled  with  8192  data  points  per  record  such  that  each  channel  of  data  consisted  of 
4096  data  points.  Finally,  the  data  were  transformed  into  desired  quantities  and  saved  as 
4-byte  floating-point  numbers.  The  header  for  this  data  file  should  consist  of  the  following 
four  numbers: 


2 32768  40  32766 


2.2.2  Remainder  of  Data  File 

Following  the  first  record  containing  the  file  header  is  the  actual  data.  When  the  data 
comes  directly  from  an  analog  to  digital  converter  (12-bit,  say),  then  the  values  range  from 
-2048  to  2047  (or  0 to  4095)  and  are  appropriately  stored  as  2-byte  integers.  Alternatively, 
if  some  preprocessing  is  typically  performed  on  the  data,  the  data  may  be  stored  as  4-byte 
floating-point  numbers.  As  a result  of  the  Fast  Fourier  Transform  (FFT)  algorithm  used  in 
this  routine,  the  number  of  data  points  per  channel  per  record  must  be  a power  of  2.  One 
or  two  separate  time  series  may  be  stored  per  record.  If  only  one  channel  of  data  is  sampled 
then  an  even  number  of  records  following  the  file  header  should  be  stored  in  the  data  file. 
This  is  necessary  because  in  the  case  of  one  channel  of  data,  two  records  are  submitted  to 
the  FFT  algorithm  to  be  transformed  simultaneously  as  a time-saving  measure.  If  an  odd 
number  of  records  is  detected,  then  the  program  automatically  appends  an  additional  record 
of  data  containing  a copy  of  the  data  values  of  the  previous  record.  When  two  channels  of 
data  are  sampled,  the  number  of  records  may  be  even  or  odd.  The  method  of  storing  the 
data  values  in  each  record  follows  the  standard  procedure  used  by  analog  to  digital  con- 
verters: 


one 


channel: 


two  channels:  Jo ^ ^2  ••• 

where  the  notation  is  as  defined  above.  This  convention  must  be  followed  even  if  the  data 
are  processed  and  then  stored  as  floating-point  numbers.  When  two  channels  of  data  are 
sampled  and  then  stored  to  a data  file,  the  values  for  each  channel  of  data  should  be  of  the 
same  data  type:  either  2-byte  integers  or  4-byte  floating-point  numbers.  If  the  user  is 
unfamiliar  with  the  task  of  creating  unformatted  FORTRAN  data  files,  then  the  user  may 
instead  create  these  files  in  an  ASCII  file  format.  A utility  MAKDAT  which  is  included  on 
the  diskette  can  then  be  used  to  write  the  file  header  and  transform  the  ASCII  data  files  into 
the  unformatted  FORTRAN  data  files  required  by  the  spectrum  routine.  See  the  section  on 
utility  files. 

2.2.3  Writing  and  reading  the  Data  File 

This  section  gives  a few  lines  of  code  which  are  representative  of  the  manner  in  which 
these  data  files  are  written  and  read.  When  writing,  the  file  is  first  opened,  the  file  header 
is  written,  the  data  is  written  and  then  the  file  is  closed. 

INTEGER*2  ICHANS,  IRSIZE,  NUMREC,  IDELTMS 

INTEGER*2  NDATA  (NMAX) 

OPEN  (2,  FILE  = ’A001.DAT’,  STATUS  = ’UNKNOWN’, 
ACCESS  = ’SEQUENTIAL’,  FORM  = ’UNFORMATTED’) 

WRITE  (2)  ICHANS,  IRSIZE,  NUMREC,  IDELTMS 


DO  J = 1,  NUMREC 
...  put  data  record  into  NDATA  array  ... 

WRITE  (2)  (NDATA  (I),  I = 1,  N) 

...  get  next  record  of  data  ... 

ENDDO 

CLOSE  (2,  STATUS  = ’KEEP’) 

The  same  procedure  is  followed  when  the  data  file  is  opened  and  read  by  the  spectrum 
routine.  The  file  is  opened,  the  file  header  is  read,  the  data  is  read  and  then  the  routine  is 
closed. 


INTEGER*2  ICHANS,  IRSIZE,  NUMREC,  IDELTMS 
INTEGER*!  NDATA  (NMAX) 

OPEN  (2,  FILE  = ’A001.DAT’,  STATUS  = ’OLD’, 

ACCESS  = ’SEQUENTIAL’,  FORM  = ’UNFORMATTED’) 

READ  (2)  ICHANS,  IRSIZE,  NUMREC,  IDELTMS 

N = IRSIZE/2 

DELT  = FLOAT  (IDELTMS)  / 1000000.0 

DO  J = 1,  NUMREC 
READ  (2)  (NDATA  (I),  I = 1,  N) 

...  process  this  record  of  data  ... 

ENDDO 

CLOSE  (2,  STATUS  = ’KEEP’) 

23  Multiple  File  Processing 

This  routine  is  capable  of  sequentially  processing  many  files  in  a batch.  The  program 
will  handle  up  to  IFMAX  files  which  are  consecutively  numbered  or  which  are  not  con- 
secutively numbered  but  have  an  arbitrary  numbering  arrangement.  All  data  files  in  the 
batch  must  start  with  the  same  first  letter.  The  letter  and  number  combination  in  the  file 
name  are  used  to  associate  each  data  file  to  a corresponding  set  of  coefficients  (if  needed) 
in  an  accompanying  coefficient  data  file;  this  is  the  reason  for  the  rather  rigid  naming 
convention.  Alternatively,  up  to  JFMAX  files  may  be  entered  on  the  command  line  when 
calling  the  program,  and  these  will  be  processed  sequentially.  These  files  need  not  have  the 
same  first  letter;  however,  these  data  files  also  cannot  use  coefficient  data  files.  This  is 
further  clarified  below.  All  data  files  must  conform  to  the  naming  convention  described 
above.  The  selection  of  options  to  perform  various  operations  on  the  data  will  then  be 
applied  to  all  files  in  the  batch. 

23.1  Consecutively  numbered  files 

If  it  is  desired  to  process  several  files  in  a consecutive  order,  then  this  option  is  selected. 
The  routine  will  ask  for  the  beginning  and  ending  file  numbers  and  then  the  first  letter  in 
the  file  names.  As  an  example,  suppose  it  was  desired  to  process  files  C031.DAT  through 
C094,DAT,  then  the  user  would  respond  as  follows: 

spectrum 

Process  files  (Qonsecutively  or  in  an  (A)rbitrary  order  : c 


Enter  starting  & ending  data  file  # : 31  94  (or  31,94) 

Enter  first  letter  of  data  file  : c 

23,2  Arbitrarily  numbered  files 

Perhaps  it  is  not  desirable  to  process  all  of  the  files  indicated  above,  but  only  selected 
ones  from  the  set.  K it  were  necessary  to  process  C031.DAT,  C042.DAT,  C067,DAT  and 
C085,DAT,  then  the  user  would  respond  as  follows. 

SPECTRUM 


Process  files  (Qonsecutively  or  in  an  (A)rbitrary  order  : A 
Enter  number  of  files  to  process  : 4 


Enter  last  digits  of  file  1:31 
Enter  last  digits  of  file  2 : 42 
Enter  last  digits  of  file  3 : 67 
Enter  last  digits  of  file  4 : 85 


Files  are  : 

31  42  67  85  (routine  echoes  file  numbers  selected) 


Enter  first  letter  of  data  file  : C 

23,3  Files  entered  on  command  line 

This  option  was  added  as  a means  of  processing  a few  files  without  having  to  answer 
all  of  the  prompts  above  and  to  allow  a bit  more  freedom  in  the  selection  of  file  names. 
Suppose  it  was  desired  to  process  TSTX.DAT,  X_IN.DAT  and  WOWl.DAT,  then  the  user 
would  enter  the  following  command  line. 


SPECTRUM  TSTX  X_IN  WOW! 


By  relaxing  the  naming  convention,  it  is  no  longer  possible  to  easily  associate  a coefficient 
data  file  with  the  correct  input  data  file.  Therefore,  if  the  user  desires  to  transform  the  integer 
data  in  each  of  the  input  files  to  voltages  only,  or  if  the  input  data  in  each  file  already  consists 
of  transformed  floating-point  numbers  (as  a result  of  some  other  program),  then  coefficient 
data  files  are  not  necessary  and  this  method  for  quickly  processing  a small  number  of  files 
may  be  used.  If  the  integer  data  in  each  input  file  must  be  transformed  into  a desired 
floating-point  quantity  using  a polynomial  transformation  prior  to  the  spectral  calculation, 


then  the  program  must  be  initiated  using  the  methods  described  in  sections  2.3.1  or  2.3.2. 
Further  discussion  on  the  manner  in  which  coefficient  data  files  are  associated  with  the 
respective  input  data  files  appears  below. 

2.4  Coefficient  Data  Files 

These  data  files  may  be  sequential  unformatted  or  sequential  formatted  FORTRAN 
data  files.  They  contain  the  constants  which  are  used  as  the  coefficients  of  the  polynomial 
that  transforms  the  input  data  from  voltages  to  some  other  desired  quantity. 

2.4.1  Naming  Convention 

The  coefficient  data  file  names  consist  of  eight  characters.  The  first  letter  must  be 
alphabetic  (A-Z)  and  must  be  the  same  letter  as  the  input  data  file  (or  series  of  input  data 
files)  with  which  it  is  associated.  The  spectrum  routine  assigns  the  remaining  seven  char- 
acters of  the  name  as  CON.DAT  (if  the  file  is  unformatted)  or  CON.PRN  (if  the  file  is 
formatted).  If  the  input  data  file  D061.DAT  were  to  be  transformed  from  voltages  into 
velocities,  then  the  coefficients  of  the  polynomial  would  be  found  in  the  coefficient  data 
file  named  DCON.DAT  (if  unformatted)  or  DCON.PRN  (if  formatted),  and  this  file  must 
reside  in  the  same  directory  as  that  of  the  executable  routine  and  the  input  data  file  unless 
specified  otherwise  in  the  configuration  file  (discussed  later).  IfD061.DAT  contained  two 
channels  of  data,  then  the  coefficients  needed  to  transform  the  data  for  the  second  channel 
would  be  found  in  the  coefficient  data  file  named  DCON2.DAT  (if  unformatted)  or 
DCON2.PRN  (if  formatted).  Again,  the  routine  requests  that  the  user  enter  the  first  letter 
only  and  the  remaining  eight  characters  are  automatically  assigned  by  the  routine. 

2.4.2  Unformatted  File  Structure 

These  coefficient  data  files  are  accessed  only  if  the  input  data  are  2-byte  integers  and 
the  user  desires  to  convert  the  integers  first  into  voltages  and  then  into  another  quantity. 
The  voltage  transformation  factor  is  not  included  in  these  files;  instead,  it  is  entered  into  the 
routine  at  the  beginning  and  therefore  applies  to  all  files  processed.  If  a series  of  files  are 
to  be  processed,  the  coefficients  for  the  transformation  from  voltages  into  the  other  quantity 
may  not  be  the  same  for  all  files  in  the  batch,  and  this  situation  is  addressed  as  follows.  The 
conversion  from  a voltage  into  a velocity,  say,  for  each  data  point  is  carried  out  using 

u„  ~ bQ  + b,v„  + b^v^  + Vn  + (2.1) 

where  v„  is  the  nth  data  point  expressed  as  a voltage,  are  the  coefficients  of 

the  velocity  transformation  and  u„  is  the  nth  data  point  which  has  been  converted  to  a 
velocity.  The  five  coefficients  of  the  transformation  remain  the  same  for  each  data  point  in 
each  record  of  the  file,  but  are  permitted  to  vary  from  file  to  file.  These  constants  are  stored 
in  a two-dimensional  array  of  4-byte  floating-point  numbers.  The  first  index  of  the  array 
is  the  number  of  the  file  to  which  the  coefficients  are  to  apply  (1  to  IFMAX),  and  the  second 


index  is  the  number  of  the  coefficient  (1=  to  5=  b^  )•  The  first  record  of  the  unformatted 
coefficient  data  file  is  a file  header  containing  two  2-byte  integer  quantities:  the  first  is  the 
number  of  sets  of  coefficients  to  read,  and  the  second  is  the  starting  file  number  to  which 
these  coefficients  should  be  associated.  The  remainder  of  the  data  file  consists  of  one  large 
record  containing  all  of  the  coefficients.  If  the  user  wishes  to  first  create  these  coefficient 
data  files  in  an  ASCII  file  format  (without  the  header),  the  utility  MAKCOEF  has  been 
included  on  the  diskette  which  will  convert  these  ASCII  files  into  the  unformatted  FOR- 
TRAN files  needed  by  the  spectrum  routine.  See  the  section  on  utility  routines.  An  example 
of  the  code  used  to  write  and  read  a sequential  unformatted  coefficient  data  file  follows. 

INTEGER*2  NUMSETS,  NUMCON,  NSTART 
REAL*4  CONST  (IFMAX,  5) 

NUMCON  = 5 

OPEN  (2,  HLE  = ’DCON.DAT’,  STATUS  = ’UNKNOWN’, 

ACCESS  = ’SEQUENTIAL’,  FORM  = ’UNFORMATTED’) 

WRITE  (2)  NUMSETS,  NSTART 

WRITE  (2)  ((CONST  (I,J),  J = 1,  NUMCON),  I = 1,  NUMSETS) 

CLOSE  (2,  STATUS  = ’KEEP’) 

The  spectrum  routine  then  reads  the  coefficient  data  file  using  code  of  the  following  form. 

INTEGER*2  NUMSETS,  NUMCON,  NSTART 
REALM  CONST  (IFMAX,  5) 

OPEN  (2,  HLE  = ’DCON.DAT’,  STATUS  = ’OLD’, 

ACCESS  = ’SEQUENTIAL’,  FORM  = ’UNFORMATTED’) 

READ  (2)  NUMSETS,  NSTART 

READ  (2)  ((CONST  (NSTART  + 1 - 1,  J),  J = 1,  NUMCON),  I = 1,  NUMSETS) 
CLOSE  (2,  STATUS  = ’KEEP’) 

The  coefficients  bQ,b^,...,b^  for  the  data  file  D061.DAT  are  then  given  by  CONST  (61,1), 
CONST  (61,2), ... , CONST  (61,5). 
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2.4.3  Formatted  Data  Files 

It  is  not  always  convenient  to  store  coefficients  in  sequential  unformatted  data  files. 
This  option  was  included  to  allow  a quick  way  of  storing  coefficients.  A sequential  formatted 
file  is  exactly  that  created  by  any  simple  text  editor.  To  create  a sequential  formatted 
coefficient  data  file  for  the  input  data  file  D061.DAT,  enter  a text  editor  or  word  processor 
(capable  of  creating  ASCII  files)  and  give  the  file  the  name  DCON.PRN.  Next,  enter  the 
five  coefficients,  •••  >^4»  one  per  line  and  then  save  the  file.  If  you  display  the  file  to 
the  screen,  it  should  look  something  like  the  following.  Note  that  the  numbers  should  all 
be  floating-point  numbers. 


3.14159265359 

(bo) 

0.0 

(bi) 

1.012345E-03 

(bi) 

0.1 

(b,) 

2.3E-05 

(b.) 

An  example  of  the  code  that  the  spectrum  routine  will  use  to  read  in  this  data  file  follows. 
REAL*4  B (5),  CONST  (IFMAX,  5) 

OPEN  (2,  RLE  = ’DCON.PRN’,  STATUS  = ’OLD’) 

READ  (2,  *)  (B  (I),  I = 1,  5) 

CLOSE  (2,  STATUS  = ’KEEP’) 


DO  I = 1,  IFMAX 
DO  J = 1,  5 
CONST  (I,  J)  = B (J) 

ENDDO 

ENDDO 

As  is  seen  from  the  above  code,  the  routine  assigns  the  five  constants  to  apply  to  all  of  the 
IFMAX  different  input  data  files  that  might  be  read.  K a different  set  of  constants  is  required 
for  each  data  file,  then  the  files  must  be  processed  one  at  a time,  or  sequential  unformatted 
data  files  as  described  above  must  be  used.  K an  input  data  file  (D061.DAT,  say)  contains 
two  channels  of  information,  then  a second  sequential  formatted  data  file  with  the  name 
DCON2.PRN  should  be  prepared  as  described  above. 
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3 Data  Preparation 

This  section  deals  with  a sequence  of  cosmetic  operations  which  may  be  applied  to 
the  data  prior  to  performing  spectral  calculations.  The  first  section  discusses  the  transfor- 
mation of  the  input  data  from  integer  values  into  voltages  and  then  optionally  into  some 
other  desired  dimensional  quantity  (for  example,  velocities).  If  the  input  data  consists  of 
floating-point  values,  then  the  transformation  operation  is  skipped.  If  desired,  the  mean 
value  may  be  removed  from  each  record  of  the  input  data,  or  a linear  trend  resulting,  perhaps, 
from  instrumentation  drift  may  be  removed  from  each  record.  Unwanted  frequencies  in  the 
input  data  can  be  removed  using  the  recursive  digital  filtering  routines  discussed  below. 
Finally,  a tapering  operation  which  may  be  coupled  with  a 50%  overlap  of  the  input  data 
records  is  included  to  eliminate  the  discontinuities  found  at  the  beginning  and  end  of  these 
records  resulting  from  the  fact  that  only  a finite  amount  of  data  can  be  sampled.  Currently, 
the  user  has  a choice  of  four  different  window  functions.  If  the  input  consists  of  two  channels 
of  data,  then  the  first  three  operations:  transformation,  mean  or  trend  removal  and  filtering 
may  be  different  for  each  of  the  two  input  channels.  The  choice  of  whether  or  not  to  employ 
tapering  as  well  as  the  choice  of  window  function,  however,  applies  to  both  input  channels. 

3.1  Transformation 

Experimental  data  sampled  by  means  of  an  analog-to-digital  (A/D)  converter  typically 
results  in  data  files  containing  an  array  of  two-byte  integer  values.  This  is  due  to  the  fact 
that  a two-byte  integer  is  capable  of  representing  numbers  in  a range  from  -32768  to  32767. 
This  range  is  sufficient  for  12,  14  and  16-bit  converters  and  is  used  to  approximate  the 
magnitude  of  an  analog  voltage  at  a given  instant  of  time.  This  voltage  is  constrained  to  lie 
within  a specified  input  voltage  range  which  may  or  may  not  be  programmable  and  which 
is  particular  to  the  A/D  converter  used.  Before  proceeding  to  the  spectral  calculation,  the 
integer  input  data  must  be  transformed  back  into  voltages  using  the  scale  factor  appropriate 
for  the  input  voltage  range  of  the  A/D  converter  used. 

During  development  of  this  routine,  a 12-bit  A/D  converter  was  used  which  was  capable 
of  digitizing  three  different  bipolar  voltage  ranges:  -IV  to  IV,  -5V  to  5V  and  -lOV  to  lOV. 
The  use  of  12  bits  allows  integer  values  ranging  from  -2048  to  2047  for  a total  of  4096 
digital  levels  to  be  used  to  quantize  each  of  the  above  input  voltage  ranges.  The  appropriate 
scale  factors  for  conversion  of  the  integer  data  back  into  voltages  are  given  by 


2V 

4096  units 


10  V 

4096  units 


and 


20  V 

4096  units 


(3.1) 


respectively.  The  user  may  choose  from  one  of  these  scale  factors  or  may  enter  an  alternative 
scale  factor  from  the  keyboard.  The  routine  first  prompts  the  user  with: 

Transform  into  (y)oltages  or  into  (O)ther  quantity 
or  read  (R)eal  data  which  is  already  transformed  : 
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If  the  input  data  consists  of  floating-point  values,  then  the  letter  R should  be  entered,  and 
the  routine  skips  the  rest  of  the  transformation  operation.  For  integer  input  data,  the  user 
responds  with  V for  transformation  into  voltages,  or  O for  voltage  transformation  followed 
by  an  additional  transformation  into  a second  dimensional  quantity;  the  following  prompt 
is  then  displayed: 

Enter  VOFST  1)  10/4096,  2)  20/4096,  3)  2/4096  or  4)  other  : 

To  this  prompt,  the  user  responds  with  a digit  (1-4)  indicating  the  desired  scale  factor  for 
voltage  transformation.  Selection  of  item  (4)  leads  to  the  prompt: 

Enter  VOFST : 

and,  at  this  point,  the  user  enters  a floating-point  number  representing  the  desired  scale 
factor. 

The  input  data  will  now  be  converted  back  into  voltages  using  the  scale  factor  specified 
above.  Since  the  analog  voltage  sampled  by  the  A/D  converter  often  comes  from  a trans- 
ducer, which  measures  some  dimensional  quantity  and  then  converts  it  to  a voltage  signal, 
provision  for  transforming  the  digitized  input  voltages  back  into  the  original  transduced 
quantity  is  included.  Knowledge  of  calibration  data  for  the  particular  transducer  is  required 
for  this  transformation,  and  the  calibration  curve  may  be  linear,  quadratic,  exponential,  etc. 
It  was  impossible  during  development  of  the  code  to  anticipate  all  forms  that  the  calibration 
curve  might  take;  therefore,  a polynomial  fit  was  decided  upon  in  hopes  that  it  would  satisfy 
the  calibration  requirements  of  most  users.  The  polynomial  may  consist  of  terms  up  to 
fourth  degree;  however,  lower  degree  polynomial  representations  may  be  used  by  setting 
the  coefficient  for  the  appropriate  term  to  zero.  The  conversion  from  a voltage  into  a velocity, 
say,  for  each  data  point  is  carried  out  using 


= ^0  + + by„ 


(3.2) 


where  foo  through  are  the  five  calibration  coefficients,  v„  is  the  nth  data  point  of  the 

current  record  which  has  just  been  transformed  into  a velocity  as  described  above  and  u„ 
is  the  resulting  nth  velocity  value.  The  coefficients  necessary  to  complete  the  transformation 
must  be  included  in  a suitable  format  in  coefficient  data  files  which  are  then  read  by  the 
routine;  instructions  for  the  creation  of  these  coefficient  data  files  are  contained  in  section 
2.4. 

An  alternative  to  the  transformation  operations  described  herein  is  also  available.  The 
user  may  transform  the  original  integer  data  from  the  A/D  converter  by  means  of  his  own 
program  into  floating-point  data  which  can  then  be  read  by  this  routine.  The  floating-point 
input  data  files  must  conform  to  the  structure  set  forth  in  section  2.  The  routine  then  skips 
the  transformations  described  above  and  proceeds  to  the  calculation  of  the  mean  and  root- 
mean-square  (rms)  values  of  each  data  record  and  averaged  values  for  the  file. 
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Summarizing  then,  the  spectrum  routine  begins  by  reading  in  an  input  data  file.  K this 
file  consists  of  integers,  then  these  integers  are  transformed  to  voltages  using  a voltage 
transformation  factor  (VOFST)  selected  by  the  user.  Optionally,  these  voltages  can  be 
further  transformed  into  some  other  dimensional  quantity.  If  the  input  data  file  contains 
floating-point  numbers,  then  the  routine  assumes  that  any  transformation  of  the  input  data 
has  already  been  applied  by  the  user.  As  the  input  data  is  read  in  by  the  spectrum  routine, 
record  by  record,  the  mean  and  rms  values  of  each  record  are  calculated  and  displayed  on 
the  screen  at  run-time.  After  all  records  of  the  input  file  have  been  read,  the  mean  and  rms 
values  for  each  record  are  averaged  to  obtain  values  which  are  representative  of  the  entire 
file.  These  averaged  mean  and  rms  quantities  are  also  displayed  on  the  screen  at  run-time. 
The  spectrum  routine  then  proceeds  to  the  various  data  preparation  operations  discussed  in 
the  next  section.  Note  that  the  averaged  mean  and  rms  values  for  the  file  are  calculated  a 
second  time  later  on  in  the  program  as  a result  of  the  spectral  computation  and  are  also 
displayed  on  the  screen.  Theoretically,  these  two  sets  of  values  should  be  identical  and 
serve  as  a check  on  the  accuracy  of  the  spectral  computation.  However,  if  any  of  the  data 
preparation  operations  to  be  discussed  next  are  applied  to  the  data  to  improve  the  accuracy 
of  the  spectral  computation,  the  two  sets  of  mean  and  rms  values  may  not  exactly  match. 
This  is  due  to  the  fact  that  the  data  preparation  operations  change  the  input  data.  The  first 
set  of  mean  and  rms  values  are  calculated  prior  to  the  data  preparation  operations,  and  the 
second  set  are  computed  from  spectral  density  calculations  which  occur  subsequent  to  the 
data  preparation  procedures. 

3.2  Mean  Removal 

For  calculation  of  spectral  quantities,  it  is  common  to  remove  the  mean  value  firom 
the  input  data.  Both  this  operation  and  the  trend  removal  operation  described  in  the  next 
section  may  be  used  for  the  removal  of  the  mean  fi’om  the  sample  data.  Trend  removal  is 
not  always  desirable;  for  that  reason  a separate  operation  to  remove  the  mean  is  needed  and 
is  described  here.  The  mean  value  for  each  record  of  data  in  the  input  file  is  calculated  firom 

X =1:1  x„  (3.3) 

n~l 

where  x„  is  defined  as  in  equation  1.1;  for  convenience,  and  to  reduce  the  algebra  in  the 

derivations  in  this  section,  n is  defined  such  that  n - 1,2,  ...,A^ , and  the  initial  point  of 
the  record  tQ  is  set  to  zero.  If  y„  is  used  to  represent  the  data  with  zero  mean,  then  this 
quantity  is  computed  from 


(3.4) 
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33  Removal  of  Linear  Trend 

As  recommended  by  Bendat  and  Piersol  (1986),  a linear  trend  removal  operation  has 
been  included  in  this  routine.  This  operation  is  designed  to  remove  low  frequency  com- 
ponents in  the  data  with  a period  longer  than  the  record  length  T . These  low  frequency 

trends  can  significantly  distort  spectral  data;  however,  trend  removal  should  only  be 
performed  if  instrumentation  drift  or  other  sources  of  trends  are  expected.  Use  of  trend 
removal  when  no  trends  are  present  can  introduce  small  contributions  to  the  spectral  density 
at  frequencies  on  the  order  of  the  reciprocal  of  the  record  length. 

This  operation  is  accomplished  by  fitting  a straight  line  through  the  data  using  the  least 
squares  fit  method.  The  values  calculated  from  this  fit  are  then  subtracted  from  the  original 
data  to  yield  data  with  zero  mean  and  with  no  frequency  components  smaller  than  l/T 
where  T is  the  record  length.  Let  ;c  be  a first  degree  polynomial  fit  to  the  original  data 
x„  (defined  as  in  equation  1.1)  where  x is  given  by 

= ^0  for  n - 1,2,...,A^  (3.5) 


and  with  and 


defined  as 


bo 


2{2N  + l)lx„ 


n - 1 


N 

6 1 nx„ 


/I  • 1 


N{N-l) 


(3.6) 


bi 


12  I nx„  - 6(A^  + 1)  i;c„ 

n “ 1 n - 1 


Ar7V'(A^-l)(A^  + l) 


(3.7) 


Let  y„  represent  the  detrended  data;  these  data  are  obtained  from  the  following  difference 


•'ft  n n 


(3.8) 


3.4  Filtering 

The  filtering  operations  incorporated  in  this  routine  were  added  to  allow  the  user  to 
pass  only  desired  frequencies  in  the  input  data  on  to  the  spectral  calculation.  Four  recursive 
(infinite  impulse  response)  Butterworth  digital  filters  for  filtering  in  the  time  domain  have 
been  included  for  this  purpose:  lowpass,  highpass,  bandpass  and  bandstop  (notch).  A short 
discussion  of  digital  filtering  in  the  time  domain  including  a summary  of  the  particular  filters 
used  in  this  routine  and  instructions  for  their  use  follows.  Note  however  that  filtering  may 
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be  applied  not  only  in  the  time  domain  but  also  in  the  frequency  domain.  When  filtering  in 
the  frequency  domain,  the  input  data  record  is  Fast  Fourier  Transformed,  then  the  FFT  output 
is  multiplied  by  a filter  function,  H{f) ; if  required,  an  inverse  FFT  could  be  performed  to 
get  the  filtered  data  back  into  the  time  domain.  Frequency  domain  filtering  may  perhaps 
be  implemented  in  a future  version  of  this  program. 

Consider  a set  of  discrete  input  values  x„  defined  as  in  equation  1.1.  Digital  filtering 

involves  a computing  algorithm  in  which  a set  of  filtered  values  y„  is  produced  as  a result 
of  operating  on  the  inputs  x„  . For  the  filters  included  in  this  routine,  this  is  a linear  operation. 
Digital  filters  can  generally  be  categorized  as  recursive  or  nonrecursive;  if  the  formula  for 
y„  contains  only  input  terms  x^  , then  the  filter  is  termed  nonrecursive.  The  general  form 
of  the  linear  nonrecursive  algorithm  is  given  by 

y.  - IKx,.,  (3.9) 

m - -M 


where  are  the  coefficients  of  the  transformation.  For  values  of  /w  < 0 , the  above 

algorithm  requires  future  values  of  the  signal  x{t) ; this  can  cause  a realizibility  problem 
in  analog  systems,  but  can  be  accomplished  in  digital  systems  where,  in  many  cases,  the 
future  values  of  x{t)  can  be  stored  in  advance.  A recursive  filter  is  one  in  which  previously 
filtered  values  such  as  y„_i , y„_2»  etc.  appear  explicitly  in  the  algorithm  for  y„ . For 
recursive  filters,  the  above  equation  is  modified  to  include  past  values  of  the  filtered  output 
y{t)  as  follows: 


M M 

y b X - y a y 

^ m n ~m  m J n 

m - -Ai  m - 1 


(3.10) 


The  transfer  function  H{f)  for  the  filter  is  related  to  the  coefficients  a„  and  b„  by 


y b„e'^^^"'^ 

H(J)  - (3.11) 

1 - 1 


This  equation  tells  how  to  determine  H(f}  from  the  coefficients  a„  and  b„  . However, 
filter  design  requires  the  inverse  problem:  how  to  obtain  a suitable  set  of  coefficients  based 
upon  knowledge  of  a desired  filter  transfer  function.  A variety  of  techniques  exist  to  solve 
this  problem,  and  the  interested  reader  is  referred  to  Steams  (1975)  for  a lengthy  discussion 
of  the  matter.  The  transfer  function  discussed  above  is  more  commonly  expressed  in  terms 
of  the  variable  z where 
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(3.12) 


For  recursive  filters,  H{z)  is  a rational  function  of  2 . It  seems  that  rational  functions  are 
particular  good  at  fitting  functions  with  sharp  edges;  most  filter  functions  are  in  this  category. 
The  filters  contained  in  this  routine  are  modified  versions  of  FORTRAN  subroutines  which 
originally  appeared  in  Steams  (1975).  The  order  of  the  filter  is  determined  by  the  number 
of  cascaded  stages  that  are  used;  a filter  of  order  P will  require  P/2  cascaded  stages.  The 
transfer  functions,  Hy^{z) , for  the  kth  stage  for  each  filter,  and  the  recursive  algorithms 
required  to  implement  those  transfer  functions  are  given  below. 


lowpass  (3.13) 


highpass  (3.14) 


A*(2'‘-2z^  + 1) 


bandpass  (3.15) 


2**  +B^Z^  + C*2^  +1^* 


“ ■^i^’/i-l  ~ Qyn-2  “ - E^yn-A 


bandstop  (3.16) 
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Note  that  the  quantity  K in  equation  3.16  above  is  a constant  which  is  the  same  for  all 
stages  of  the  notch  filter.  The  other  coefficients,  through  E,, , are  a more  compact 
notation  for  the  a„  and  b„  values  for  each  stage  k . In  general,  through  change 
from  stage  to  stage  and  are  evaluated  from  complicated  expressions  which  the  interested 
reader  may  find  in  the  source  code. 

To  use  these  filters,  the  user  selects  the  desired  filter  type  and  enters  the  critical  (-3  dB) 
frequency  for  the  low  and  highpass  filters.  For  the  bandpass  and  bandstop  filters,  two  critical 
frequencies  must  be  specified.  The  number  of  filter  sections  to  cascade  must  also  be  spe- 
cified. The  order  of  the  filter  is  twice  the  number  of  cascaded  sections.  As  the  order  of  the 
filter  is  increased,  the  transfer  function  for  the  filter  becomes  sharper  and  more  closely 
approximates  the  ideal;  however,  the  amount  of  computation  time  is  also  significantly 
increased.  A good  starting  point  may  be  to  use  five  cascaded  sections  (tenth  order);  this  is 
a reasonable  trade-off  between  accuracy  and  computational  efficiency.  These  filters  work 
well  when  the  input  data  is  not  deterministic;  however,  unlike  nonrecursive  filters,  recursive 
filters  can  become  unstable  and  produce  an  output  which  grows  exponentially.  This  has 
been  observed  in  a few  rare  instances  during  testing  of  this  program  and,  so  far,  has  been 
found  to  occur  only  when  the  input  data  has  been  a perfect  computer-generated  sinusoid. 
The  user  should  be  aware  of  this  limitation. 
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3.5  Tapering 

The  tapering  option  has  been  added  to  the  routine  as  a grooming  operation  in  order  to 
improve  the  resulting  spectral  density  estimates.  A wide  variety  of  window  functions  may 
be  found  in  the  literature  including  extended  discussions  concerning  their  performance 
characteristics.  See,  for  example,  Bendat  and  Piersol  (1986)  and  Press,  etal.  (1986).  Rather 
than  repeat  these  long  discussions  here,  a brief  summary  of  the  need  for  tapering  will  be 
followed  by  a description  of  the  window  functions  available  for  use  with  this  routine. 

As  discussed  earlier  in  regard  to  equation  1.1,  the  data  processed  by  this  routine  is  in 
the  form  of  a finite  set  of  samples  x„  of  the  continuous  signal  x{t)  which  exists  for  a time 
T =N  At . Now,  suppose  that  x{t)  is  a finite  portion  of  a continuous  signal  v(r)  which  is 
of  unlimited  extent  in  time.  Then,  x{t)  can  be  considered  to  be  the  product  of  v(t)  and  a 
rectangular  window  function  w(f) 

x(0  = w(r)v(r)  (3.17) 


where 


|1  0 

[O  otherwise 


(3.18) 


A plot  of  the  rectangular  (square)  window  function  w(f)  is  shown  in  figure  3.1  where,  for 
convenience,  T is  chosen  to  equal  one.  The  computation  of  spectral  density  estimates  for 
x{t)  involves  the  calculation  of  the  Fourier  transform  of  x(t)  which  is  equivalent  to  the 
Fourier  transform  of  the  product  of  w{t)  and  v{t) . A theorem  from  Fourier  analysis  relates 
the  Fourier  transform  of  a product  of  two  functions  in  time  to  the  convolution  of  the 
transforms  of  these  functions.  In  other  words,  the  following  is  a Fourier  transform  pair: 

w{t)v{t)  W{f)  * V{f)  (3.19) 

where 

W(f)*V(f)  . H/(a)  V{f  - a)  da  (3.20) 


The  spectral  density  of  interest  would  be  exact  if  the  transform  V(f)  could  be  calculated. 
Instead,  X(f)  = W(f)  * V{f)  is  the  quantity  which  is,  in  fact,  calculated.  This  problem  then 
is  seen  to  be  directly  related  to  the  finite  length  of  the  data  to  be  transformed  and  to  the  fact 
that  the  data  is  sharply  truncated  at  each  end  of  the  time  record.  A plot  of  the  modulus  of 
the  Fourier  transform  of  the  rectangular  window  function,  W{f) , may  be  found  in  figure  3.2 
and  is  defined  as 


|W(f)/W(0)  I 
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Figure  3.1  Window  functions  available  for  tapering  the  data 


XT  XT 

f 

Figure  3.2  Normalized  amplitude  of  Fourier  transform  of  window  functions 
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(3.21) 


Again,  for  plotting  convenience,  T is  chosen  to  equal  one.  This  function  is  characterized 
by  a large  main  lobe  with  smaller  side  lobes  which  slowly  decrease  in  size.  The  ideal  window 
function  would  be  a delta  function;  that  is,  an  infinitely  thin  main  lobe  with  no  side  lobes 
at  all.  As  W{f)  is  convolved  with  V{f) , information  at  frequencies  well  separated  from 
the  main  lobe  can  leak  into  the  calculation  and  significantly  distort  the  resulting  spectra. 
The  large  amplitude  sidelobes  need  to  be  reduced  to  suppress  this  problem. 

This  is  accomplished  by  multiplying  the  original  data  by  a window  function  which 
tapers  the  time  history  at  the  beginning  and  end  of  the  record  to  remove  the  discontinuities 
there.  The  tapering  functions  included  in  this  routine  are  the  most  commonly  used  ones, 
and  the  choice  of  a particular  function  requires  a trade-off  between  the  thickness  of  the  main 
lobe  (resolving  power)  and  the  size  of  the  side  lobes.  The  following  window  functions  (the 
subscripts  merely  indicate  the  name  of  the  window)  may  be  selected: 


otherwise 


Hanning 


(3.22) 
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(3.23) 


0 


otherwise 


\ 2 
0 


0:sr  :sr 
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(3.24) 


otherwise 
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= 


1 - 

2 

T 

2 

0 

O^t^T 

otherwise 


Parzen 


(3.25) 


These  functions  are  plotted  in  figure  3.1.  The  moduli  of  the  Fourier  transforms  for  each  of 
the  above  windows  are  plotted  in  figure  3.2  and  are  normalized  by  their  value  at  / = 0 . The 
fact  that  these  transforms  do  not  equal  one  for  /-O  implies  that  the  resulting  spectral 
estimates  would  be  reduced  in  magnitude  upon  application  of  one  of  these  windows.  To 
avoid  this  problem,  the  appropriate  scale  factor  is  computed  by  the  routine  when  any  of 
these  windows  are  selected  and  automatically  applied  to  the  data  before  computing  the 
spectral  estimate  in  order  to  preserve  the  correct  amplitudes.  The  Fourier  transform  for  each 
of  the  above  windows  is  given  below. 
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In  equations  3.26  and  3.27  above  = 1/7’.  Recall  that  ifnoneofthese  windows  are  selected 
(no  tapering),  then  the  square  window  defined  in  equation  3.18  and  with  a Fourier  transform 
as  defined  in  equation  3.21  is  implicitly  used.  The  most  commonly  used  window  function 
is  the  Hanning  window;  the  transform  for  this  function  has  the  smallest  side  lobes  but  also 
has  the  widest  main  lobe.  Also  shown  is  a window  which  provides  a cosine  taper  for  the 
first  and  last  tenths  of  the  data  record,  is  one  for  the  remainder  of  the  data  record  and  is 
denoted  by  the  name  Tukey  (also  known  as  a 10%  cosine  taper).  This  window  has  been 
used  by  many  researchers  in  the  hopes  that  it  will  throw  away  less  of  the  data.  Examination 
of  the  modulus  of  the  Fourier  transform  for  this  window  shows  that  it  manages  to  narrow 
the  main  lobe  a small  amount  relative  to  the  other  windows,  but  it  does  this  at  the  expense 
of  large  side  lobes  that  are  nearly  as  large  as  that  of  the  square  window  (no  tapering  at  all). 

A better  way  to  throw  away  less  of  the  data  is  to  use  any  of  the  windows  available 
above  and  then  to  use  overlapped  processing  techniques.  Before  describing  overlapped 
processing,  it  is  useful  to  review  the  procedure  for  the  calculation  of  spectra  when  no  overlap 
is  used.  The  data  are  prepared  (detrended,  filtered,  etc.),  and  then  the  data  are  tapered.  The 
data  are  tapered  record  by  record  and  transformed  to  get  the  necessary  spectra.  When  50% 
overlapped  processing  is  used,  the  tapering  operation  proceeds  a bit  differently.  First,  the 
initial  record  of  the  data  is  tapered.  Then,  the  second  half  of  the  first  record  is  grouped 
together  with  the  first  half  of  the  second  record  to  form  a new  second  record,  and  then  this 
new  record  of  data  is  tapered.  Then,  the  two  halves  of  the  old  second  record  form  the  new 
third  record,  and  this  is  then  tapered.  Then,  the  second  half  of  old  second  record  is  grouped 
together  with  the  first  half  of  the  old  third  record  to  form  a new  fourth  record  which  is  then 
tapered.  This  proceeds  throughout  the  file  until  M old  records  are  processed  creating 
exactly  2M  new  records.  The  last  grouping  pairs  the  second  half  of  the  old  last  record  with 
the  first  half  of  the  old  first  record  to  form  the  new  2Mth  record.  This  last  record  was 
included  in  order  to  simplify  the  code  for  this  procedure.  Because  the  2Mth  record  has  a 
discontinuity  at  the  midpoint,  this  record  should  not  be  transformed.  However,  simply 
removing  this  record  would  require  that  an  odd  number  of  records  be  submitted  to  the  FFT 
code.  As  pointed  out  earlier,  the  FFT  routine  has  been  designed  such  that  it  transforms  two 
records  of  data  in  one  pass  as  a time-saving  measure;  this  requires  (in  the  one  channel  case) 
that  an  even  number  of  records  be  processed  by  the  FFT  routine.  Therefore,  the  overlapping 
operation  has  not  been  modified;  2M  records  are  created  from  M original  records  for  both 
the  one  and  two  channel  cases.  These  2M  records  are  then  submitted  to  the  FFT  subroutine 
and  are  transformed;  but  the  results  of  the  2Mth  record  are  discarded  after  the  transformation. 
Summarizing,  when  overlapping  is  employed,  2M  records  are  created  and  are  then  trans- 
formed, but  only  the  results  from  2M  - 1 records  are  actually  used  in  the  final  calculations. 
To  avoid  discontinuities  resulting  from  the  use  of  overlapping  in  any  of  the  other  2M  - 1 
records,  the  M original  records  must  be  contiguous  in  time.  The  data  need  not  be  contiguous 
from  record  to  record  if  overlapping  is  not  employed.  The  use  of  50%  overlapped  processing 
doubles  the  number  of  tapering  and  FFT  operations;  however,  this  will  retrieve  about  90% 
of  the  stability  lost  due  to  the  tapering  operation  and  will  decrease  the  variability  of  the 
resulting  estimate. 
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Finally,  tapering  is  best  used  when  the  input  data  are  essentially  random  data  — that 
is,  not  deterministic.  When  this  program  was  tested,  input  data  (both  deterministic  and 
random)  with  known  spectral  densities  were  used.  When  a sinewave  was  used  for  input, 
the  power  and  amplitude  spectra  exhibited  a lower  noise  floor  when  no  tapering  was  used. 
On  the  other  hand,  when  colored  noise  (bandwidth  limited  white  noise)  was  used  for  input, 
a substantial  improvement  was  evident  in  the  output  spectra  when  any  of  the  four  window 
functions  described  above  was  used  to  taper  the  data. 
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4 Spectrum  Computation 

This  section  introduces  the  definitions  of  the  discrete  Fourier  transform  and  the  various 
spectral  density  functions  computed  by  this  routine.  The  means  by  which  these  estimates 
are  computed  as  well  as  the  errors  associated  with  the  calculation  are  discussed.  The  need 
for  additional  code  to  adjust  the  value  of  the  computed  phase  angle  of  the  cross-spectral 
density  function  is  explained  by  means  of  an  example,  and  an  explanation  of  the  coherence 
function  calculated  by  this  routine  is  given. 

4.1  Fast  Fourier  Transform 

Using  the  notation  of  Bendat  and  Piersol  (1986),  the  Fourier  transform  of  a function 
x{t) , where  x:(0  may  be  real  or  complex-valued,  is  given  by 


X{f)  = 


c 


dt 


J 


.00 


(4.1) 


This  definition  is  not  particularly  useful  for  experimentally  derived  data  which  exist  for  only 
a finite  time  interval.  For  this  reason,  one  usually  employs  a finite-range  transform  over 
the  interval  O^t  , and  this  transform  is  defined  by 


(4.2) 


As  in  equation  1.1,  assume  that  each  record  of  the  input  data  file  consists  of  N samples  of 
x{t)  which  have  been  obtained  at  equally  spaced  points  in  time  with  spacing  At . Note  that 
to  has  been  set  equal  to  zero  for  convenience.  Then,  each  record  contains  the  values  x„ 
where 


x„  = x(nAt)  for  /i  = 0, 1,2,  ...,iV  - 1 


(4.3) 


Now,  the  finite-range  transform  defined  in  equation  4.2  can  be  written  in  discrete  form  as 


X{f,T)  = At 


N-l 


1 


(4.4) 


Since  only  N input  values  of  x„  are  used  in  the  computation  of  the  transform,  only  N 
output  values  of  X(J,T)  may  be  specified.  These  are  typically  chosen  at  the  following 
discrete  frequency  values 
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T NAt 


(4.5) 


Thus,  the  discrete  form  of  the  finite-range  Fourier  transform  becomes 

AT-l 

X(J,J)^At  2 x„  e =AtX^  for  A: -0,l,2,...,Ar-l  (4.6) 

rt  » 0 

Although  N output  values  have  been  specified,  results  are  unique  only  out  to  k =-N  12  ; 
this  is  a result  of  the  symmetry  of  the  Fourier  transform  for  real  input  values.  Thus,  the 
highest  frequency  that  can  be  resolved  for  a sampling  rate  of  1/Ar  samples  per  second  is 
the  Nyquist  frequency  given  by 


4 = /c  = 


1 

2 At 


(4.7) 


which  occurs  for  k = N 12  . To  compute  the  X^  values  in  equation  4.6  requires  complex 
multiply-add  operations  for  each  record  of  the  data  file;  note  that  1 complex  multiply-add 
operation  is  equivalent  to  4 real  multiply-adds.  This  can  be  quite  a formidable  computational 
task  as  N becomes  large.  Fortunately,  several  algorithms  for  reducing  the  complexity  and 
increasing  the  speed  of  the  computation  exist  and  are  known  collectively  as  fast  Fourier 
transforms  (FFT).  Rather  than  discuss  these  algorithms  in  detail,  the  reader  is  referred  to 
the  discussions  that  may  be  found  in  the  references.  The  FFT  algorithm  used  here  employs 
a technique  known  as  time  decomposition  with  input  bit  reversal.  The  kernel  of  this  algorithm 
was  taken  from  Steams  (1975).  The  simplest  implementation  of  this  algorithm  calls  for  N 
to  be  a power  of  2.  This  is  required  in  this  routine,  and  if  necessary,  each  record  of  the  input 
data  file  should  be  padded  with  zeros  at  the  end  to  satisfy  this  requirement. 

The  FFT  algorithm  used  in  this  routine  is  designed  to  handle  the  general  case  where 
x{t)  is  a complex-valued  function;  that  is,  the  input  data  which  is  passed  to  the  FFT  subroutine 
consists  of  two  arrays:  one  containing  the  real  parts  and  the  other  containing  the  imaginary 
parts.  However,  experimentally  sampled  random  time  series  data  are  real-valued.  When 
the  FFT  of  a purely  real-valued  function  is  computed,  certain  symmetries  in  the  resulting 
transform  allow  one  to  compute  the  FFTs  for  two  real-valued  functions  simultaneously  by 
inserting  one  into  the  input  array  of  real  parts  and  the  other  into  the  array  of  imaginary  parts. 
The  trick  is  to  unscramble  the  individual  transforms  from  the  resulting  computation.  This 
speed  improvement  has  been  incorporated  into  this  routine,  but  it  places  a restriction  on  the 
input  data  files.  If  the  input  data  file  to  be  processed  by  SPECTRUM  consists  of  one  channel 
of  sampled  data,  then  the  data  file  must  consist  of  an  even  number  of  records.  The  reason 
for  this  is  that  input  data  records  are  submitted  to  the  FFT  algorithm  two  at  a time  to  obtain 
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the  speed  improvement  discussed  above.  If  the  input  data  file  consists  of  two  channels  of 
sampled  data,  then  each  record  of  the  file  is  first  decomposed  into  the  two  separate  channels, 
and  these  are  then  submitted  to  the  FFT  algorithm  to  be  simultaneously  processed.  So,  for 
the  case  of  two  channels  of  sampled  data,  the  input  data  file  may  have  an  even  or  odd  number 
of  records. 

4.2  Autospectral  (Power)  Density 

The  one-sided  autospectral  density  function,  defined  for  positive  frequencies  only,  is 
usually  expressed  as  twice  the  amplitude  of  the  Fourier  transform  of  the  autocorrelation 
function  as  follows: 


GJf)  - 2r^RJj)e-^f'dz 

for  0 ^ ^ oo  , otherwise  zero  . 


(4.8) 


where  RJ^t)  is  the  autocorrelation  function  defined  as 


RJj)  = F[x(04^+x)]  (4.9) 

and  the  notation  E[  ] denotes  the  expected  value  operation.  GJ(f)  is  a real  function  of  / 
regardless  of  whether  x{t)  is  real  or  complex.  This  calculation  is  implemented  in  this 
routine  using  digital  computing  techniques  employing  finite  fast  Fourier  transforms  of  the 
original  data  records.  Using  this  method,  the  autospectral  density  may  be  obtained  from 


GJJ)  - 2 lim  ^ E[\X(f,T)\^]  . 

T — oo  1 


(4.10) 


I I denotes  the  modulus  of  a complex  quantity,  and  X{f,T)  is  the  finite  Fourier  transform 
which  is  valid  for  the  interval  O^t  and  is  calculated  as  in  equation  4.6.  If  the  input 
data  file  contains  two  channels  of  data,  denoted  as  x{t)  and  y{t) , then  selection  of  the 
power  spectral  density  option  results  in  the  computation  of  G^(f)  and  Gyyif) . 

Specifically,  the  calculation  proceeds  by  computing  the  FFT  for  each  record  of  the 
input  data  file,  determining  the  modulus  of  the  complex  result,  squaring  it  and  then  saving 
the  results  to  an  array.  This  is  repeated  for  each  record;  upon  completion,  the  results  are 
ensemble  averaged  to  approximate  the  expected  value  operation  required  in  the  definition. 
The  computation  of  an  expected  value  is  needed  to  reduce  the  random  error  of  the  estimated 
autospectral  density  function.  The  number  of  records  in  the  input  data  file  determines  the 
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number  of  members  in  the  ensemble  which  will  be  averaged.  If  the  number  of  records  in 
the  data  file  (before  tapering)  is  denoted  by  , then  the  normalized  random  error  of  the 
autospectral  density  estimate  is  given  by 


(4.11) 


where  the  normalized  random  error  is  defined  as  a firactional  quantity  by  dividing  the  standard 
deviation  of  the  estimate  by  the  actual  value  estimated.  If  the  tapering  operation  is  used,  it 
will  increase  the  variability  of  the  resulting  estimate  by  an  amount  which  depends  upon  the 
particular  window  function  chosen.  This  increase  in  the  variance  of  the  estimate  is,  in  the 
worst  case  (Hanning),  about  two;  but  if  the  50%  overlap  technique  is  also  employed,  it  will 
counteract  this  increase  and  will  recover  about  90%  of  the  variability  caused  by  tapering. 
Thus,  the  normalized  random  error  calculated  firom  equation  4.11  will  serve  as  a reasonably 
accurate  estimate  for  both  the  no  taper  and  taper  with  overlap  cases.  If  tapering  is  used  but 
overlapping  is  not,  then  the  variance  is  doubled  and  the  normalized  random  error  defined 
in  equation  4.11  should  be  multiplied  by  the  square  root  of  two.  The  routine  computes  the 
appropriate  value  of  the  normalized  random  error  and  displays  it  on  the  screen  at  run-time. 
As  an  example,  consider  an  autospectral  density  estimate  computed  from  a data  file  with 
100  original  records  which  are  then  overlapped  and  tapered  to  produce  199  records;  this 
estimate  will  have  a normalized  random  error  of  10%.  To  reduce  this  error  to  5%,  the  input 
data  file  would  have  to  contain  400  records.  There  is  no  specific  numerical  limit  to  the 
number  of  records  that  an  input  data  file  may  contain;  the  spectrum  routine  has  successfully 
computed  spectral  densities  for  input  data  files  with  10,000  records.  However,  data  files 
with  more  than  100  records  are  usually  quite  large.  Furthermore,  the  spectrum  routine  must 
create  temporary  files  during  computation,  and  the  size  of  these  files  (depends  upon  the  size 
of  the  input  data  file)  is  often  very  large.  Reserving  a sufficient  amount  of  space  on  the  hard 
disk  or  on  a RAM  disk  for  these  files  is  typically  the  limiting  factor  on  the  number  of  records. 

The  computation  of  a fast  Fourier  transform  and  power  spectral  density  function  allows 
the  calculation  of  the  mean  and  root-mean-square  (rms)  values  of  each  record  of  the  input 
time  series.  The  equation  for  the  mean  value  of  each  data  record  can  be  obtained  by  setting 
/:  = 0 in  equation  4.6.  One  then  obtains 


N Nn.o 


The  mean-square  value  of  each  record  can  be  computed  from 


(4.12) 


'^GJf)df 


= E[x\t)]  = 


0 


(4.13) 
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The  rms  value  is  just  the  square  root  of  this  number.  The  spectrum  routine  calculates  the 
mean  and  rms  values  of  the  input  time  series  using  this  procedure  and  reports  them  on  the 
screen  at  run-time  along  with  the  appropriate  value  of  the  normalized  random  error. 

43  Cross-spectral  Density 

4.3.1  Calculation  Procedure 

A one-sided  cross-spectral  density  function  (f)  may  be  defined  in  a manner  similar 
to  equation  4.8;  this  quantity  is,  in  general,  complex  and  is  given  by 

G^(/)  = 2jy4x)e-‘^^dr  - C^(f)  - i Q^{J)  (4.14) 

for  0 ^ / :£  00  , Otherwise  zero  . 


The  quantity  is  the  cross  correlation  function  defined  by 

Rxyi'^)  = E[xit)y(t +x)]  (4.15) 

The  cross-spectral  density  may  also  be  computed  using  fast  Fourier  transform  procedures 
and  is  the  method  adopted  here.  The  equivalent  expression  to  equation  4.10  for  the  one-sided 
cross-spectral  density  function  is 

G^(f)  - 2 lim  i E[r(f,T)  Y(f,T)]  (4.16) 

7-.  00  i 


where  the  notation  X*{f,T)  is  used  to  represent  the  complex  conjugate.  Alternatively,  this 
can  be  expressed  in  polar  notation  as 


G^(f) 

where 


0 S/<  00 


G:^(f)  \ = 
e^CO  = 


K(f)  + Qlif)] 
1 Q^if) 


tan 


C^(f) 


(4.17) 


and 
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Note  that  this  option  may  be  selected  only  when  the  input  data  file  contains  two  channels 
of  data,  then,  the  magnitude  and  phase,  | (f)  \ and  0^  (f)  respectively,  are  the  quantities 

which  are  calculated  and  then  saved  in  the  output  file. 

As  described  above  in  the  section  on  power  spectral  density  calculation,  it  is  possible 
to  compute  the  value  of  the  autocorrelation  function  for  x =-  0 (the  mean-square  value)  by 
calculating  the  integral  of  the  autospectral  density  function.  In  a similar  manner  one  may 
compute  the  cross  correlation  function  evaluated  at  x » 0 from 

(‘‘•18) 

The  spectrum  routine  calculates  this  quantity  and  displays  it  on  the  screen  at  run-time. 
4.3.2  Adjustment  of  Phase  Angle 

In  order  to  clarify  the  need  for  additional  code  to  adjust  the  values  of  the  phase  angle 
^xyif)  > example  taken  from  Hess  (1990)  will  be  given  which  also  serves  to  illustrate  the 
use  of  the  cross  spectral  density  function  in  general.  Specifically,  a point  measurement  of 
the  vertical  displacement  of  a compliant  surface,  responding  to  the  unsteady  forcing  of  a 
fluid  flow  above  the  surface,  was  acquired;  call  this  displacement  signal  x(t) . As  x(t)  was 
measured,  a second  displacement  signal  was  simultaneously  measured;  call  this  second 
signal  y(t) . The  locations  on  the  compliant  surface  of  the  measurement  of  x(t)  and  y(t) 
were  as  follows:  both  were  obtained  on  the  centerline  of  the  rectangular  compliant  slab, 
and  y(r)  was  located  11.4  cm  downstream  of  x(t) . The  fluid  flow  was  turbulent;  the 
instantaneous  streamwise  component  of  velocity  above  the  compliant  surface  consisted  of 
a small  fluctuating  quantity  superimposed  upon  a larger  mean  value.  A disturbance  in  the 
fluid  flow  caused  a displacement  of  the  compliant  surface  which  was  measured  at  the 
upstream  location  and  was  recorded  at  time  t as  x{t) . As  the  flow  disturbance  was 
convected  downstream  by  the  mean  flow,  the  displacement  of  the  compliant  surface  also 
travelled  downstream  and  was  measured  at  the  downstream  location  as  y{t  +x)  where  x 
was  the  time  delay  for  the  disturbance  to  travel  11.4  cm.  Quantities  of  interest,  then,  are 
the  convection  speed  of  the  displacement  fluctuation  on  the  compliant  surface  as  well  as 
the  decrease  in  amplitude  of  this  fluctuation  as  it  decays. 

For  this  particular  situation,  the  decay  of  the  displacement  field  may  be  modeled  by 
the  time  delay  problem,  see  Bendat  and  Piersol  (1986).  This  model  provides  a method  for 
determining  the  convection  speed  of  the  displacement  fluctuation.  Consider  two  zero-mean 
stationary  random  signals  x{t)  and  y(t)  where  y(t)  is  defined  by 


y{t)  = ax(t  -Xq)  + n{t)  . 


(4.19) 


The  value  a is  a constant  attenuation,  Tq  is  a constant  time  delay  and  n{t)  is  uncorrelated, 
zero  mean  value  noise  at  the  output.  This  process  is  illustrated  in  figure  4.1. 


n(t) 


y(t) 


The  cross-correlation  of  these  two  signals  can  be  computed  to  be 


(4.20) 


Thus,  the  cross-correlation  is  an  attenuated  replica  of  the  autocorrelation,  which  is  displaced 
in  time  by  the  amount  Tq  . Without  going  through  all  of  the  details,  a one-sided  cross-spectral 
density  function  (f)  may  also  be  computed;  this  quantity  is,  in  general,  complex  and 
for  this  problem  is  given  by 


(4.21) 

with  \G^{f)\  = aG^if)  and  Q^(f)  = ImJ  . 

Since  the  time  delay  Tq  appears  only  in  the  phase  angle,  0^(/) , it  can  be  determined  from 

measurement  of  the  phase  angle  and  by  noting  that  it  is  a linear  function  of  the  fi-equency 
with  slope  2icZq  . Summarizing  then,  a cross-spectral  density  is  computed  for  x(t)  and 
y{t) . If  the  process  can  be  modelled  by  the  time  delay  problem,  then  the  modulus  will  be 
an  attenuated  replica  of  the  autospectral  density,  G^(f) , and  the  phase  angle  will  be  a linear 
function  of  the  frequency  with  a slope  given  by  2kZq  . A convection  speed  may  then  be 

11  4 

calculated  from  11^=  — cm/sec. 

An  example  of  this  procedure  is  shown  in  figure  4.2.  The  top  graph  shows  the  auto- 
spectral density  curve  G^{f)  and  below  it  the  modulus  | G^(f)  | for  the  corresponding 
cross-correlation;  indeed,  the  bottom  curve  does  appear  to  be  an  attenuated  copy  of  the  first. 
Note  the  substantial  additional  scatter  in  the  cross-spectral  density  estimate  relative  to  that 
for  the  autospectral  density  estimate  at  the  higher  frequencies.  The  normalized  random  error 


32 


f(H2) 

Figure  4.2  Cross-spectral  density  modulus  example 

for  the  cross-spectral  density  estimate  turns  out  to  be  a function  of  frequency  unlike  that  for 
the  autospectral  density  estimate.  A discussion  of  the  means  for  estimating  this  error  is 
given  in  a later  section. 

Figure  4.3  displays  the  computed  phase  angle  expressed  in  radians.  With  the  data 
presented  in  this  form,  it  is  difficult  to  fit  the  data  with  a least  squares  fit  in  order  to  determine 
the  slope,  and  therefore,  Tq  . This  difficulty  may  be  traced  to  the  fact  that  an  inverse  tangent 
must  be  computed  in  order  to  calculate  the  value  of  the  phase  angle.  The  FORTRAN  intrinsic 
function  ATAN2  returns  a value  for  the  inverse  tangent,  but  forces  the  result  to  lie  in  the 
range  -n  s result  ^ n . However,  the  phase  angle  function  is  periodic  with  period  2:i . 
Using  this  fact,  additional  code  may  be  used  to  straighten  out  the  curve  by  adding  integer 
multiples  of  2k  as  required.  Figure  4.4  shows  the  adjusted  values  of  the  phase  angle  and 
a least  squares  fit  to  the  slope.  The  calculated  values  of  Tq  and  / U* , where  is  the 
mean  flow  speed  far  above  the  compliant  surface,  are  given  on  the  plot.  The  scatter  in  the 
phase  angle  data  becomes  large  at  a frequency  for  which  the  energy  in  the  signal  is  about 
two  decades  below  the  peak  value.  Computed  values  of  the  convection  speed  using  this 
method  vary  by  less  than  2%  from  those  obtained  using  an  alternative  method. 

The  phase  angle  data  can  be  adjusted  manually,  but  this  is  a tedious  and  time-consuming 
procedure.  Instead,  the  spectrum  routine  automatically  computes  both  the  original  values 
of  the  phase  angle  as  well  as  the  adjusted  values  of  the  phase  angle  when  the  cross  spectral 
density  and  the  coherence  function  are  calculated.  Note  that  the  relationships  derived  in 
equations  4.21  apply  only  when  equation  4. 19  may  be  used  as  a model  for  the  data.  Therefore, 
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both  sets  of  phase  angle  data,  unadjusted  and  adjusted  values,  are  stored  to  the  output  file 
and  made  available  to  the  user.  The  task  of  adjusting  these  phase  angle  values,  while 
conceptually  simple,  is  somewhat  tricky  to  program.  A description  of  the  technique  follows. 

The  standard  deviation  of  the  phase  angle  estimate  as  a function  of  frequency  is 
computed  when  the  coherence  function  is  computed  (this  will  be  explained  in  the  section 
on  error  estimates).  The  average  of  a moving  window  of  NSD  values  of  this  standard 
deviation  is  computed  and  checked  against  a threshold  value  given  by  SDLIM  . If  the 
average  standard  deviation  is  too  high,  then  the  phase  angle  points  are  considered  to  be 
essentially  random  and  no  further  attempts  to  adjust  the  values  are  made.  Thus,  when  this 
condition  occurs,  the  remaining  phase  angle  values  are  kept  the  same  as  the  original  values 
(not  adjusted).  When  the  moving  average  of  the  standard  deviation  is  less  than  SDLIM  , 
then  adjustment  may  take  place  as  follows.  A straight  line  is  fit  to  a moving  window  of 
NPTS  values  of  the  original  phase  angle  data  using  a least  squares  fit.  The  coefficients  of 
this  fit  are  used  to  guess  the  next  value  of  the  phase  angle;  if  the  actual  value  of  the  next 
phase  angle  data  point  differs  from  the  guessed  value  by  an  amount  greater  than  n radians, 
then  multiples  of  2k  radians  are  added  to  (or  subtracted  from)  the  actual  value  until  the 
difference  between  the  actual  and  guessed  values  is  less  than  k radians.  The  procedure 
continues  to  adjust  subsequent  values  of  the  phase  angle  until  the  moving  average  of  the 
standard  deviation  exceeds  SDLIM  . This  portion  of  the  code  contains  three  adjustable 
parameters:  NSD  = 7 , NPTS  = 7 and  SDLIM  = 70.0®  . These  values  have  been  selected 
by  trial  and  error  to  yield  the  best  response  for  the  data  most  often  processed  by  the  author. 
Other  values  may  be  chosen  by  changing  the  numbers  in  the  appropriate  PARAMETER 
statement  at  the  beginning  of  the  code.  Note  that  NSD  and  NPTS  must  be  chosen  such 
that 


NSD  ^ 2 NPTS  - 1 (4.22) 

is  satisfied.  This  adjustment  procedure  seems  to  work  reasonably  well  for  most  situations 
encountered  by  the  author;  however,  some  manual  adjustment  of  values  may  still  be  required 
in  some  circumstances. 
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Figure  4.3  Phase  angle  example 


Figure  4.4  Phase  angle  - adjusted  values 


43.3  Coherence  Function 

The  coherence  function,  also  known  as  the  coherency  squared  Junction  is  defined  by 
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I G^U)  1^ 

GM  G„{f) 


(4.23) 


where  0 ^ ^%{f)  ^ 1 for  all  / . This  function  is  analogous  to  the  square  of  the  more  well 
known  correlation  coefficient  function  given  by 


, (t-)  ^ E[xit)y{t +t)] 


(4.24) 


where  -1  s P^(t)  ^ 1 for  all  values  of  x . As  described  by  Bendat  and  Piersol  (1986),  the 

function  P;^(x)  measures  the  degree  of  linear  dependence  between  x{t)  and  y{t)  for  a 
displacement  of  x in  y{t)  relative  to  x{t) . Similarly,  for  linear  systems,  the  function 
ylyif)  measures  the  fractional  portion  of  the  mean  square  value  at  the  output  y{t)  that  is 
contributed  by  the  input  x{t)  at  frequency  /.  On  the  other  hand,  the  quantity  [1  -Y^(/)] 
is  a measure  of  the  mean  square  value  of  y{t)  not  accounted  for  by  x(t)  at  frequency  f . 
Ifjc(r)  and  _y(r)  are  completely  unrelated,  the  coherence  function  will  be  zero.  Furthermore, 
Bendat  and  Piersol  (1986)  indicate  that  if  the  coherence  function  is  greater  than  zero  but 
less  than  unity,  one  or  more  of  the  following  three  possible  physical  situations  exist. 

1.  Extraneous  noise  is  present  in  the  measurements. 

2.  The  system  relating  x(r)  and  y{t)  is  not  linear. 

3.  y{t)  is  an  output  due  to  an  input  jc(r)  as  well  as  to  other  inputs. 

The  calculation  of  this  function  requires  not  only  the  computation  of  the  cross  spectral 
density  function  but  also  the  computation  of  power  spectral  densities  for  each  channel  of 
data.  The  code  has  been  optimized  to  minimize  the  amount  of  additional  computational 
effort;  however,  in  anticipation  of  the  possibility  that  a slow  computer  might  be  employed, 
the  routine  prompts  the  user  for  this  option  at  run-time. 

43.4  Error  Estimates 

The  definition  of  the  coherence  function  makes  possible  at  this  point  a discussion  of 
the  error  associated  with  estimates  of  the  magnitude  and  phase  of  the  cross  spectral  density 
function  and  of  the  coherence  function  itself.  This  information  comes  from  Bendat  and 
Piersol  (1986).  The  normalized  random  error  of  the  estimate  of  the  magnitude  of  the  cross 
spectral  density  function  is  given  by 
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O^if) 


(4.25) 


where  | Yj^(/)  | is  the  positive  square  root  of  the  coherence  function  evaluated  at  / , and 

is  the  number  of  records  in  the  input  data  file.  The  error  approaches  that  for  the  power 
spectral  density  estimate  (see  equation  4.1 1)  as  | y^(f)  | goes  to  one,  and  the  error  is  higher 
for  values  less  than  one.  Unlike  the  error  estimate  for  the  power  spectral  density  function, 
all  of  the  error  estimates  in  this  section  are  a function  of  frequency. 

Since  the  phase  angle  may  assume  the  value  zero,  a normalized  random  error  cannot 
be  specified;  instead  the  standard  deviation  (s.d.)  of  the  phase  angle  estimate  is  given  by 


s.d.if) 


I Ixyif)  I 


(4.26) 


The  standard  deviation  is  seen  to  approach  zero  as  | Y^(/)  | goes  to  one  and  does  so  inde- 
pendent of  for  n^>\  . Continuing,  the  normalized  random  error  of  the  coherence 
function  itself  can  also  be  calculated  and  is  defined  as 


E(/)  = 


I YxyC/)  I 


(4.27) 


where  the  estimated  value  of  the  coherence  function  is  substituted  as  an  approximation  to 
the  unknown  true  values  of  Y^(/)  ^nd  | Y;o-(/)  I called  for  in  equation  4.27.  Note  that  the 
normalized  random  error  approaches  zero  as  either  oo  or  y^  — »•  1 . This  implies  that 
coherence  function  estimates  can  be  more  accurate  than  either  autospectra  or  cross  spectra 
estimates  (which  are  used  to  calculate  the  coherence  function)  when  Y^(/)  is  close  to  unity. 
If  the  coherence  function  is  calculated,  the  spectrum  routine  automatically  computes  the 
three  error  estimates  defined  in  equations  4.25  to  4.27  as  a function  of  frequency  and  saves 
the  results  in  a separate  output  file.  The  file  name  and  structure  of  this  file  will  be  discussed 
in  the  section  on  output  files.  Equations  4.25  to  4.27  are  valid  for  the  no  taper  and  taper 
with  overlap  cases.  K tapering  is  used  without  overlapping,  then  the  error  estimates  in 
equations  4.25  to  4.27  must  each  be  increased  by  multiplying  by  the  square  root  of  two,  and 
this  is  automatically  performed  by  the  spectrum  routine. 

It  is  perhaps  useful  to  know  the  number  of  records  required  in  the  input  data  file  to 
attain  a specified  normalized  random  error  or  standard  deviation  for  each  of  the  estimates 
mentioned  above.  This  has  been  calculated  for  various  values  of  Y^(/)  for  the  cases: 
^ [ I |]  = 0.1 , E[Yi(/)]=0.1  and  5.  i/.  [ 0^(/)]  = 0.052  - 3°  , and  the  results 

are  given  in  table  4.1  below. 


37 


0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

e[|G.,(r)|]-0.1 

334 

250 

200 

167 

143 

125 

112 

426 

274 

183 

122 

79 

46 

21 

e[yI(/)]=0.1 

327 

180 

100 

54 

26 

10 

2 

Table  4.1  Number  of  records  for  specified  error 
4.4  Amplitude  Spectral  Density 

The  amplitude  spectrum  option  was  added  to  allow  the  user  to  examine  the  Fourier 
coefficients  directly.  The  amplitude  spectrum  is  defined  as 

F,(/)  - lim  £:[|XC/-,r)|]  (4.28) 

7 00 

where  | X{f,  T)  \ is  the  modulus  of  the  quantity  calculated  in  equation  4.6.  If  the  input  data 
file  contains  two  channels  of  data,  denoted  as  x{t)  and  y{t) , then  selection  of  this  option 
results  in  the  computation  of  FJJ)  and  F^if) . 

A simple  example  taken  from  Steams  (1975)  is  used  to  demonstrate  the  utility  of  this 
option.  Consider  the  following  function  which  was  computer-generated  and  then  saved  as 
a set  of  digital  values  to  simulate  the  sampling  process: 


x(r)  = e~*  sin  t (4-29) 

A plot  of  a portion  of  this  function  is  shown  in  figure  4.5.  The  Fourier  transform  of  this 
function  can  be  derived  using  the  definition  in  equation  4.1  as 


X{f)  = 


1 

(/  2jt/  +1)^  + 1 


and  the  modulus  of  this  quantity  becomes 


(4.30) 


1 


\m\  - 


(4.31) 


la)  di 


Figure  4.6  Plot  of  amplitude  spectrum  calculation 
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This  input  function  is  deterministic,  and  therefore  only  one  record  of  data  would  normally 
be  required;  however,  because  an  even  number  of  records  are  required  by  the  spectrum 
routine  for  one  channel  of  input  data,  two  identical  records  of  this  function  were  generated. 
N = 2048  samples  per  record  were  created  with  a time  interval  of  At  = 0.01  seconds  between 
successive  samples.  A portion  of  the  results  is  shown  in  figure  4.6  with  the  symbols  rep- 
resenting the  values  computed  by  the  spectrum  routine  and  the  continuous  line  the  values 
generated  by  equation  4.31. 
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5 Output  Data  Files 

For  every  input  data  file  processed,  the  spectrum  routine  creates  one  (or  sometimes 
two)  ASCn  output  files  containing  the  results,  A second  output  file  is  created  only  if  the 
user  processes  two  channels  of  input  data  and  computes  the  cross-spectral  density  and  the 
coherence  function;  the  second  output  file  then  contains  error  estimates  of  cross-spectral 
magnitude  and  phase  and  of  the  coherence  function. 

5.1  Naming  Convention 

The  output  file  name  will  consist  of  nine  characters.  The  last  four  characters  specify 
the  file  extension;  the  default  file  extension,  [.PRN],  will  be  used  unless  specified  otherwise 
in  the  configuration  file.  The  first  four  characters  of  the  output  file  name  will  match  the 
first  four  characters  of  the  corresponding  input  data  file.  The  fifth  character  is  used  to 
designate  the  type  of  file  that  is  being  written.  The  four  possibilities  are: 


power  spectrum  P 

cross  spectrum  C 

amplitude  spectrum  A 

error  estimates  (if  computed)  E 


The  output  file  will  be  written  to  the  same  location  as  the  input  data  file;  this  will,  by  default, 
be  the  location  where  the  executable  version  of  the  spectrum  routine  resides.  This  location 
can  be  changed,  however,  in  the  configuration  file. 

5.2  File  Structure 

These  output  data  files  consist  of  ASCII  data  as  opposed  to  binary  data.  They  may  be 
easily  edited  by  any  ASCII  editor,  imported  into  spreadsheet  programs  for  further  manip- 
ulation or  imported  into  plotting  programs  to  be  graphed.  These  files  contain  numbers  only; 
no  labels  are  included  to  describe  the  various  columns  of  data  in  the  output.  While  somewhat 
inconvenient  for  the  user,  this  was  done  to  prevent  the  confusion  that  sometimes  occurs 
when  data  containing  both  numbers  and  labels  are  imported  into  other  software  programs. 
In  the  sections  that  follow,  a thorough  description  of  the  output  to  be  expected  for  each  of 
the  four  types  of  output  files  is  given,  and  this  should  allow  the  user  to  interpret  the  results 
without  any  difficulty. 

5.2.1  Power  Spectrum  Output  Files 

The  first  line  of  the  output  contains  the  mean  and  rms  values  computed  as  described 
in  equations  4.12  and  4.13  and  then  ensemble-averaged  over  all  records  of  the  file.  Note 
that  the  mean  and  rms  values  reported  here  are  calculated  subsequent  to  any  cosmetic 
operations  that  may  have  been  applied  to  the  data  to  improve  the  quality  of  the  spectral 
density  computation.  Removal  of  the  mean  value,  detrending,  filtering  and  tapering  oper- 
ations may  modify  the  mean  and  rms  so  that  they  do  not  match  the  original  values  for  the 
file  which  were  calculated  and  reported  at  run-time  prior  to  the  application  of  these 
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operations.  Continuing,  the  next  line  of  the  output  file  is  blank.  The  calculated  results  begin 
on  the  third  line  and  consist  of  four  or  six  columns  of  data  depending  on  whether  the  input 
file  contains  one  or  two  channels.  The  quantities  are  stored  in  the  following  order; 

One  channel: 

mean  rms 

skip  a line 

/ log/  G^if)  \og  G^if) 

1 

NAt 

1 

2At 

Two  channels: 

mean  rms 

skip  a line 

/ log/  GM  ^ogG^if)  G^(f)  log  G^(f) 

1 

NAt 

1 

2At 


In  either  case  the  output  file  contains  N /2  rows  of  data  after  the  blank  line  with  each  row 
calculated  for  values  of  / ranging  from  I /NAt  to  l/2Ar  with  a resolution  of  A/=  1/NAt  . 
This  structure  is  typical  for  all  of  the  output  files  created  by  the  spectrum  routine,  and  in  the 
following  descriptions  of  the  output  only  the  column  labels  will  be  identified. 

5.2.2  Cross  Spectrum  Output  Files 

A cross  spectrum  may  only  be  calculated  when  the  input  file  contains  two  channels 
of  data.  However,  the  output  file  is  somewhat  different  depending  upon  whether  or  not  the 
coherence  function  is  calculated.  In  either  case,  the  file  begins  with  the  value  of  the  cross 


correlation  function  evaluated  at  x = 0 and  calculated  according  to  equation  4.18.  The  next 
line  is  skipped,  and  the  following  record  is  the  first  oi  N H rows  of  data  arranged  in  five 
or  seven  columns  as  described  below. 
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Coherence  function  not  calculated: 

skip  a line 

/ log/  o^(/)  logG^cn  e^(f) 

Coherence  function  calculated: 

^^(0) 
skip  a line 

/ log/  G^if)  logG^C/)  e^(/)  ylif) 


5.2.3  Amplitude  Spectrum  Output  Files 

When  the  amplitude  spectrum  option  is  chosen,  the  structure  of  the  output  file  is  slightly 
different  than  for  the  other  files.  No  mean  or  rms  value  is  reported.  Instead,  the  file  contains 
N /2  + 1 rows  of  data  arranged  in  the  columns  shown  below.  The  first  row  consists  of  data 
computed  for  / = 0 ; data  for  this  value  of  / is  not  stored  when  other  options  are  selected, 
but  is  useful  when  examining  the  Fourier  coefficients.  Actually,  the  value  / = £ is  used; 
this  is  necessary  .since  log(O)  is  undefined.  The  value  of  £ is  defined  in  a PARAMETER 
statement  at  the  beginning  of  the  code  and  is  currently  set  at  £ = 1.0  :c  10'^  . 

One  channel: 

/ log/  F,C0 

£ 

1 

NAt 


2At 
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Two  channels: 

/ log/  F^(f)  F^if) 

E 

1 

NAt 


1 

2Ar 


5.2.4  Error  Estimate  Output  Files 

When  both  the  cross-spectral  density  and  the  coherence  function  options  are  selected, 
the  spectrum  routine  calculates  the  errors  associated  with  these  estimates  as  a function  of 
frequency  using  equations  4,25  through  4.27.  These  error  data  are  then  stored  in  a separate 
file  from  the  other  results.  This  file  contains  N /2  rows  of  data  arranged  in  the  columns 
identified  below.  The  columns  containing  normalized  random  errors  consist  of  dimen- 
sionless numbers,  and  the  standard  deviation  of  the  phase  angle  is  expressed  in  degrees. 


/ log/  £[|g,(/)|]  s.d.[e.,(/)]  e[y^(/)] 
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6 Run-time  Considerations 

Two  topics  are  discussed  in  this  section  which  assume  importance  when  the  program 
is  executed.  The  creation  of  a configuration  file  allows  the  user  to  bend  some  of  the  rules 
regarding  file  name  creation  and  the  location  of  data  files  as  well  as  temporary  files  created 
by  the  spectrum  routine.  The  second  topic  concerns  the  manner  in  which  the  size  of  the 
temporary  files  created  by  the  routine  may  be  determined. 

6.1  Configuration  File 

The  presence  of  an  ASCII  file  with  the  name  SPECTRUM. CFG  in  the  same  directory 
as  the  SPECTRUM.EXE  file  allows  optional  settings  to  take  effect  which  may  have  a 
favorable  impact  on  the  execution  speed  and  overall  ease  of  use  of  the  routine.  This  ASCII 
file  allows  the  user  to  specify  three  parameters  which  will  become  default  values  when  the 
program  is  executed. 

1.  The  first  item  is  the  drive  and  path  information  for  the  location  of  input,  output 
and  coefficient  data  files.  If  the  user  wishes  to  specify  that  these  files  may  be  found 
in  (and  written  to)  the  same  location  as  the  SPECTRUM.EXE  file,  then  the  word 
ROOT  should  be  entered  on  a line  by  itself  Otherwise,  the  desired  drive  and  path 
information  should  be  entered  in  the  following  format:  C:\SPECTRUM\DATA\ 
Note  that  the  last  backslash  is  required. 

2.  The  second  item  in  the  file  is  the  drive  and  path  information  which  specifies  where 
temporary  data  files  created  by  the  spectrum  routine  will  be  stored.  If  a RAM  disk 
is  present  on  the  user’s  computer  and  is  large  enough  to  hold  the  temporary  files, 
then  the  execution  speed  of  the  routine  will  be  greatly  enhanced  if  these  files  are 
located  on  this  RAM  disk.  A means  for  estimating  the  size  of  these  temporary  files 
is  discussed  below.  Again,  the  user  specifies  ROOT  or  the  desired  drive  and  path 
information  using  a format  such  as:  D:\TEMP\ 

3.  The  last  item  allows  the  user  to  change  the  file  extension  for  output  data  files  from 
its  default  value  of  [.PRN]  to  some  other  desired  value.  The  user  should  place  the 
extension  on  a line  by  itself  in  the  following  format:  .OUT 

Note  that  the  leading  period  is  required. 

Comment  lines  may  be  included  in  the  configuration  file;  these  lines  must  start  with  an 
asterisk  character.  The  routine  is  not  case-sensitive  --  the  options  may  be  specified  in  upper 
or  lower  case.  The  spectrum  routine  does  not  require  the  presence  of  the  configuration  file; 
if  this  file  is  absent,  then  the  following  default  values  are  used:  1.  root  2.  root  3.  [.pm]  If 
a configuration  file  is  used,  then  it  must  reside  in  the  same  directory  as  the  SPECTRUM.EXE 
file  and  all  three  options  must  be  specified.  An  example  configuration  file  follows. 
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* This  is  a sample  configuration  file  for  Spectrum.  Comment  lines 

* and  blank  lines  must  begin  with  an  asterisk. 

* 

* The  first  item  in  the  file  is  the  path  where  the  input  data 

* files  are  found  and  to  which  the  output  data  will  be  written. 

* Acceptable  values  are:  root,  c:\ , g:\ , c:\test\bin\ , etc.  The 

* path  must  end  with  a backslash  if  root  is  not  used, 
root 

* The  second  item  in  the  file  is  the  drive  on  which  the  temporary 

* files  are  to  be  created.  Program  execution  speed  will  increase 

* dramatically  if  this  is  a RAM  disk. 

* Acceptable  values  are:  root,  c:\ , g:\ , c:\test\bin\ , etc.  The 

* path  must  end  with  a backslash  if  root  is  not  used. 
d:\ 

* The  third  item  in  the  file  is  the  default  file  extension  to  use 

* for  the  output  data  files.  If  not  specified  in  this  file,  then 

* the  .PRN  extension  will  be  used  in  order  that  the  data  may  be 

* imported  into  LOTUS  for  plotting.  The  extension  must  consist  of 

* a period  followed  by  three  letters. 

.prn 

6.2  Temporary  Files 

The  input  data  files  that  are  to  be  processed  by  the  spectrum  routine  are  typically  very 
large;  therefore,  the  data  is  analyzed  one  record  at  a time.  Furthermore,  the  final  result  of 
the  calculation  is  not  arrived  at  in  one  pass;  instead,  intermediate  results  are  generated  and 
stored  in  temporary  files  which  are  then  used  to  obtain  the  final  results.  These  temporary 
files  are  automatically  deleted  upon  completion;  however,  sufficient  space  must  be  reserved 
on  the  hard  disk  or  on  a RAM  disk  for  temporary  use  by  these  files.  The  following  discussion 
explains  how  these  files  are  created  and  allows  the  user  to  estimate  their  size  based  upon 
the  size  of  the  input  data  file.  Note  that  if  several  input  files  are  to  be  processed  by  the 
routine,  only  the  temporary  files  for  the  input  data  file  currently  being  processed  reside  on 
the  hard  disk  at  any  one  time;  the  user  needs  to  allow  sufficient  space  based  upon  the  size 
of  the  largest  input  data  file. 

Two  temporary  files  are  created  for  each  channel  of  input  data.  After  a record  of  the 
input  data  file  is  read  in  by  the  routine,  the  integer  data  is  sorted  into  one  or  two  arrays 
depending  on  the  number  of  input  channels;  each  channel  is  then  transformed  into  floating 
point  data  (if  necessary),  the  mean  value  or  a linear  trend  is  removed  and  any  desired  filtering 
takes  place.  These  data  are  then  stored  to  temporary  files  with  the  extensions  : [.TRF]  - 
channel  1 and  [.TR2]  - channel  2 (if  required).  This  process  is  repeated  until  all  records  of 
the  input  data  file  are  read.  If  the  input  data  file  consisted  of  2-byte  integers,  then  the  [.TRF] 
file  will  be  twice  the  size  of  the  input  [.DAT]  file  since  the  [.TRF]  file  contains  4-byte 
floating  point  numbers.  If  the  input  data  file  consists  of  two  channels  of  integer  data,  then 
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the  [.TRF]  and  [.TR2]  files  will  each  be  the  same  size  as  the  [.DAT]  file  (the  [.TRF]  and 
[.TR2]  files  each  contain  half  as  many  numbers  as  the  input  [.DAT]  file,  but  each  number 
is  a 4-byte  value).  If  the  input  data  file  consists  of  one  channel  of  4-byte  real  numbers,  the 
[.TRF]  file  will  be  the  same  size  as  the  [.DAT]  file,  and  finally,  for  two  channels  of  4-byte 
values  in  the  input  file,  the  [.TRF]  and  [.TR2]  files  will  each  be  half  the  size  of  the  [.DAT] 
file.  After  creating  the  [.TRF]  and  [.TR2]  (if  required)  files,  the  spectrum  routine  proceeds 
to  the  tapering  operation.  If  tapering  is  selected,  the  routine  reads  a record  of  data  at  a time 
from  the  [.TRF]  and  [.TR2]  files,  creates  two  new  records  for  each  channel  (if  overlapping 
is  used)  and  stores  the  results  in  two  temporary  files  with  the  extensions  [.TAP]  and  [.TP2]. 
The  last  step  in  the  process  is  to  transform  the  data  stored  in  the  [.TAP]  and  [.TP2]  files  and 
perform  the  required  spectral  calculation.  The  final  data  is  written  to  one  output  file 
regardless  of  the  number  of  channels  of  input  data  and  the  extension  [.PRN],  or  alternatively 
an  extension  specified  in  the  configuration  file  is  used.  K tapering  is  not  selected,  then  the 
[.TAP]  and  [.TP2]  files  are  not  created,  and  storage  space  does  not  need  to  be  allocated.  If 
the  [.TAP]  and  [.TP2]  files  are  created,  each  will  be  twice  the  size  of  the  previous  [.TRF] 
or  [.TR2]  file,  since  two  tapered  records  are  written  for  every  input  record  read  to  implement 
the  50%  overlap  technique.  (If  tapering  is  employed  but  overlapping  is  not,  then  the  [.TAP] 
and  [.TP2]  files  will  be  the  same  size  as  the  [.TRF]  or  [.TR2]  file.  The  table  belows  gives 
the  worst  case  scenario  and  assumes  overlapping  is  used  when  tapering  is  used.)  At  the  end 
of  the  calculation,  only  the  original  [.DAT]  file  and  the  final  [.PRN]  file  will  remain.  If  no 
alternative  is  specified  in  the  configuration  file,  the  temporary  files  will  be  created  in  the 
same  directory  as  the  input  data  files.  Several  examples  follow  which  show  the  sizes  of  the 
various  files  based  upon  an  input  data  file  size  of  1 Mb. 


Two  Channels  ? 

No 

Yes 

No 

Yes 

No 

Integers  ? 

Yes 

Yes 

No 

No 

Yes 

Tapering  ? 

Yes 

Yes 

Yes 

Yes 

No 

[.DAT] 

1 Mb 

1 Mb 

1Mb 

1Mb 

1Mb 

[.TRF] 

2 Mb 

1 Mb 

1 Mb 

500  kb 

2Mb 

[.TR2] 

— 

1 Mb 

— 

500  kb 

— 

[.TAP] 

4 Mb 

2Mb 

2Mb 

1 Mb 

— 

[.TP2] 

— 

2 Mb 

— 

1Mb 

— 

[.PRN] 

64  kb 

93  kb 

48  kb 

83  kb 

64  kb 

Total  Temp 
space  needed 

6Mb 

6Mb 

3Mb 

3Mb 

2Mb 

Table  4.2  Memory  required  for  temporary  files 
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7 Examples 

Three  simple  examples  are  considered  that  serve  to  illustrate  various  features  of  the 
spectrum  routine.  The  example  using  colored  noise  demonstrates  the  use  of  tapering  and 
allows  the  calculated  power  spectrum  to  be  compared  with  a known  power  spectrum.  The 
computation  of  the  amplitude  spectrum  of  a sinewave  reveals  the  importance  that  sampling 
parameters  can  have  on  the  computed  result,  and  the  last  example  shows  the  effects  of  the 
application  of  a variety  of  digital  filters  prior  to  the  computation  of  power  spectra.  A few 
comments  regarding  the  considerations  necessary  for  the  acquisition  of  experimental  data 
files  are  given  at  the  end. 

7.1  Colored  Noise 


A routine  to  produce  bandwidth-limited  white  noise  (colored  noise)  has  been  written 
and  used  extensively  during  the  testing  phase  of  the  development  of  the  spectrum  routine. 
This  program  prompts  the  user  for  a few  parameters,  then  generates  the  colored  noise 
sequence  and  stores  the  output  data  in  a manner  that  simulates  the  sampling  process.  This 
data  may  then  be  input  into  the  spectrum  program,  and  the  output  power  spectrum  estimate 
may  be  compared  with  the  known  power  spectrum  for  colored  noise.  This  handy  routine 
is  a useful  utility  program  and  has  been  included  on  the  diskette  with  the  name  COLOR. 

The  procedure  for  producing  colored  noise  begins  by  generating  random  numbers  that 
vary  from  0 to  1.  This  sequence  is  uniformly  distributed  over  the  interval  from  0 to  1.  Let 
x„  be  the  nth  random  number,  then  x„  - 1/2  will  be  uniformly  distributed  over  the  interval 
from  -1/2  to  1/2  . Now,  recall  that  the  mean  and  variance  of  a uniformly  distributed 
random  variable  is  given  by 


a + b 

and 


(b-af 

12 


(7.1) 


where  a and  b are  the  endpoints  of  the  interval.  Using  these  definitions,  we  find  that  the 
mean  and  variance  of  the  sequence  x„- 1/2  is 


= 0 and 


(7.2) 


If  a power  spectrum  for  this  sequence  were  calculated,  then  one  would  obtain  the  mean 
value,  , the  mean-square  value,  ip;  and  the  function  G^{f) . Note  that  = cr,  + p;  , 
so  that  for  this  particular  sequence 


G^(f)  df 


total  power  (P) 


(7.3) 
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If  these  data  are  considered  to  be  samples  of  a continuous  signal  x(t)  with  sampling  interval 
A/  , then  the  power  spectrum  should  look  like  that  shown  in  figure  7.1. 


Gxxtn 


2PAt 


Figure  7.1  White  noise  spectrum  with  total  power  P 


For  convenience,  it  is  useful  to  adjust  the  total  power  in  the  sequence  in  order  that  the  power 
density  G^(/)  = 1 for  0 where  = l/2Ar  is  the  Nyquist  frequency.  To  accomplish 

this,  consider  the  sequence 


(7.4) 


where  x„  is  a random  number  on  the  interval  from  0 to  1 as  before  and  y„  ranges  from 
-K 12  to  KH  . For  this  sequence,  the  mean  and  variance  are 


p,  = 0 


and 


o"  = — 
" 12 


(7.5) 


Now,  choose  P = l/2Ar,  so  that 


of  - GJf)  df 


12  2At 


(7.6) 


From  this,  we  find  that  -K”  = , and  for  this  value  of  K , the  power  density  for  the 

sequence  defined  in  equation  7.4  satisfies 


1 

0 otherwise 


(7.7) 
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A plot  of  this  power  spectrum  is  shown  in  figure  7.2. 


G„(f) 


Figure  7.2  White  noise  spectrum  with  total  power  l/2Af 

The  final  step  is  to  create  a new  sequence  by  passing  the  values  y„  through  a bandpass 
filter  with  transfer  function  Hif) . The  power  density  for  this  new  sequence  is  given  by 


GJJ)  GJf) 


(7.8) 


A recursive  digital  bandpass  filter  with  \H{f)  | = 1 and  described  in  equation  3.15  has  been 
used  in  the  COLOR  routine.  The  user  specifies  the  two  (-3  dB)  points,  fy  and  /2 , and  the 
number  of  cascaded  filter  stages  to  use.  The  order  of  the  filter  is  twice  the  number  of  cascaded 
stages.  The  output  of  the  program  is  a sequence  z„  with  a power  spectrum  as  shown  in 
figure  7.3. 


Figure  7.3  Colored  noise  spectrum  with  total  power  (/2  -/i) 
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The  color  routine  was  used  to  create  the  sequence  z„  (simulated  sampling  of  the 

continuous  signal  z{t) ) which  was  stored  in  a data  file  that  could  be  read  by  the  spectrum 
routine.  The  file  contained  100  records  with  2048  data  points  per  record.  The  spacing 
between  data  points  was  0.01  seconds.  The  parameters  for  the  bandpass  filter  were: 
/i  = 20  Hz  , /2  = 30  Hz  and  five  cascaded  stages  for  a tenth  order  filter.  A plot  of  a portion 
of  the  simulated  time  series  appears  in  figure  7.4.  The  data  file  was  processed  by  the  spectrum 
routine  several  times  using  various  choices  of  tapering  functions,  and  the  power  spectra  are 
plotted  in  figure  7.5.  When  computing  these  spectra,  no  mean  removal  or  detrending  was 
used,  and  no  further  filtering  was  performed.  As  can  be  seen,  GJJ)  - 1 in  the  passband 
and  rolls  off  to  zero  very  rapidly  elsewhere.  (Note  that  the  base  10  logarithm  of  these 
numbers  are  plotted  in  order  to  magnify  the  rolloff  region.)  This  plot  clearly  shows  the 
advantages  gained  with  the  use  of  tapering  and  overlapped  processing  as  opposed  to  no 
tapering;  however,  in  this  case,  the  choice  of  tapering  function  makes  little  difference. 
Table  7.1  shows  the  rms  values  reported  for  each  case;  all  of  these  are  within  0.6  % of  the 
theoretical  value,  and  the  minor  variations  are  due  to  the  fact  that  the  rolloff  is  not  perfect. 


0 0.2  0.4  0.6  0.8  1 

t (sec) 

Figure  7.4  Simulated  time  series  z{t) 


Theory 

No  Taper 

Hanning 

Parzen 

Tukey 

Welch 

rms  values 

3.16228 

3.14463 

3.14665 

3.14617 

3.14315 

3.14429 

error  (%) 

— 

0.56 

0.49 

0.51 

0.60 

0.57 

Table  7.1  RMS  values  reported  by  spectrum  routine 


Figure  7.5  Power  spectra  for  colored  noise;  total  power  = VlO 


f(H2) 

Figure  7.6  Effect  of  detrend  on  power  spectra 
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For  data  generated  by  the  computer,  no  linear  trend  in  the  data  is  to  be  expected,  and  the 
detrending  operation  should  not  be  used.  The  indiscriminate  use  of  this  operation  can  lead 
to  added  power  at  / = 1 ///Ar  since  this  operation  is  applied  record  by  record;  this  can  only 
occur  if  the  power  at  this  frequency  is  lower  than  the  small  amount  of  noise  generated  by 
the  use  of  this  operation.  In  almost  every  case  where  the  random  data  originates  from  an 
experiment  (vis-a-vis  computer  generated  perfect  random  data)  this  problem  will  not  occur, 
but  the  user  should  be  aware  of  the  possibility.  A plot  of  the  added  power  that  occurs  with 
the  use  of  the  detrending  operation  is  shown  in  figure  7.6. 

The  following  is  an  example  session  showing  the  various  prompts  and  responses  that 
were  used  with  the  color  routine  for  this  example. 


color 

Enter  the  low  frequency  cutoff  : 20.0 
Enter  the  high  frequency  cutoff : 30.0 
Enter  delta  T seconds.  (1.0/Sample  Rate) : 0.01 
Enter  # of  filter  sections  to  cascade  : 5 
Enter  # of  points  per  record  (N)  : 2048 
Enter  # of  records  (NUMREC)  : 100 
Enter  a negative  integer  : -5 
Enter  output  file  name  (4  chars) : COLR 


(Defines  passband) 

(Sampling  rate  = 100  Hz) 

(Tenth  order  filter) 

(Must  be  power  of  two) 

(No  limit  on  # records) 

(Required  seed  value  for) 
(random  number  generator) 


COLORED  NOISE  UTILITY 


Creating  COLR.DAT  now. 

Writing  record  0100 

Use  voltage  transformation  factor  = 20.0  / 4096.0 
Program  terminated  successfully. 


Next  is  an  example  session  showing  the  prompts  and  responses  needed  to  process  the 
COLR.DAT  file  using  the  spectrum  routine. 


spectrum  coir 

Transform  into  (V)oltages  or  into  (O)ther  quantity 
or  read  (R)eal  data  which  is  already  transformed  : v 

Enter  VOFST  1)  10/4096,  2)  20/4096,  3)  2/4096  or  4)  other  : 2 

Remove  (M)ean  value,  (D)etrend  or  (N)either  : n 

Filter  the  data  using  (L)owpass,  (H)ighpass 

(B)andpass,  Band(S)top  or  (N)one  : n 

How  many  channels  of  data  (1  or  2)  : 1 

(A)mplitude,  (P)ower  or  (Qross  spectrum  : p 

Taper  the  data  (Y/N)  : y 

(H)anning,  (P)arzen,  (T)ukey  or  (W)elch  window  : h 
Use  overlap  (Y/N)  : y 

Processing  COLR.DAT  now. 

Record  0100  Ave  = .15974E-03  Rms  = 3.1035 

File  Totals  : Ave  = .26703E-05  Rms  = 3.1446 

Taper  applied  = Hanning. 

50%  overlap  used. 

0200  records  will  be  written. 

Tapering  record  0200. 

Transforming  records  0199  and  0200. 

Record  0200  was  discarded. 

Writing  COLRP.PRN  now. 

mean  rms 

2.18332E-06  3.1467 

norm,  random  error  = .10000 

starting  time  : 09:15:01.71 

ending  time  : 09:16:06.68 

elapsed  time  : 00:01:04.97 


Program  terminated  successfully. 
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Note  that  the  maximum  number  of  points  per  record  is  controlled  by  the  value  of 
NMAX.  The  value  of  NMAX  is  currently  set  to  16384,  (If  the  user  should  ever  desire  to 
make  a change  to  the  value  of  NMAX,  the  change  must  be  made  both  in  the  COLOR  routine 
and  also  in  the  SPECTRUM  routine.  NMAX  must  be  the  same  number  in  both  routines  at 
all  times.) 

7.2  Sinusoidal  Input 

When  developing  the  spectrum  routine,  it  was  felt  that  a simple  test  would  be  to  attempt 
to  determine  the  Fourier  coefficient  (amplitude)  of  a sinewave  by  computing  its  amplitude 
spectrum.  Although  this  can  be  done,  it  is  not  a trivial  procedure  and  is  therefore  repeated 
here  as  an  example.  A sinewave  generation  routine  was  written  and  is  included  on  the 
diskette  as  SINE;  this  routine  creates  a data  file  that  can  be  read  by  the  spectrum  routine. 
Each  record  of  this  data  file  contains  a portion  of  the  time  series 


x{t,k)  = A sin[2ji^f  + 9*] 


(7.9) 


where  A is  the  amplitude,  is  the  frequency  and  0*  is  a phase  angle  which  changes  from 
record  to  record  and  is  a function  therefore  of  the  record  number  k . The  presence  of  a 
changing  phase  angle  from  record  to  record  is  merely  to  simulate  sampling  of  noncontiguous 
chunks  of  the  continuous  function.  This  phase  angle  will  not  alter  the  value  of  the  amplitude 
spectrum  and  is  therefore  left  out  of  the  following  derivation  for  simplicity. 

The  finite-range  Fourier  transform  defined  in  equation  4.2  can  be  computed  for  this 
function  x{t)  and  is  given  by 


X{fJ)  = 


2i 


sin  Ji  (/ - /o)  T -iK{f-f^)T 


sinjt(/  + /o)r  -iK(f*f,)T 
n(f+fo)T  ^ 


(7.10) 


This  function  must  now  be  evaluated  at  / = /o , so  taking  the  limit  gives 


AT 

lim  X{f,T)  = — 

/-*/o 


sin2jl/or  -iTnf^T 

“ 2nfoT  ^ 


(7.11) 


where  the  first  term  was  evaluated  using  THopital’s  rule.  The  modulus  of  this  quantity  is 
a function  of  the  record  length  T , and  is  given  by 


AT 

lim  X{f,T)  \ - — 


sin^  2ji  /o  r sin  AnfaT 


{2KfoTf 


2ji/o  T 


+ 1 


(7.12) 


To  further  examine  this  dependence  upon  the  record  length,  one  can  express  T in  terms  of 
n , the  number  of  cycles  of  the  sinewave,  using 
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T = NAt 

fo 


(7.13) 


With  this  value  for  the  record  length,  equation  7,12  becomes 


lim  X{f,T)\  = 


2 


sin^  2nji 
{2n7if 


sin  Ann 
2nn 


(7.14) 


The  function  /(n)  is  plotted  in  figure  7.7  and  is  shown  to  oscillate  about  the  value  one. 
The  extent  of  the  oscillation  decreases  with  increasing  n , and  in  the  limit  as  the  record 
length  goes  to  infinity,  the  function  equals  one. 


n 


Figure  7.7  Oscillatory  function  f{n) 

Due  to  the  finite  nature  of  the  signal,  one  can  obtain  the  limiting  value  for  the  Fourier 
coefficient  only  by  paying  careful  attention  to  the  choice  of  record  length.  Specifically,  if 
the  record  length  is  chosen  to  contain  an  integral  number  of  cycles  of  the  sinewave,  then  as 
can  be  seen  from  figure  7.7,  f{n)  = l for  integer  « . Summarizing  then,  the  limiting  value 
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of  the  amplitude  spectrum  of  the  sinewave  can  be  computed  if  the  record  length  is  chosen 
to  contain  an  integral  number  of  cycles  of  the  sinewave;  when  this  criterion  is  satisfied,  then 
the  modulus  is  given  by 


lim  X(J,T)\  = — for  integer  n 

/-/o  2 


(7.15) 


Now,  using  the  definition  given  in  equation  4.28,  the  value  of  the  amplitude  spectrum 
reported  by  the  spectrum  routine  at  the  frequency  /q  will  be  the  ensemble  average  over  all 
records  of  the  data  file  of  the  modulus  defined  in  equation  7.15.  Three  further  constraints 
are  placed  upon  the  record  length  in  addition  to  the  requirement  of  integral  values  of  n : N 
must  be  a power  of  two,  N should  be  large  for  better  accuracy  and  At  which  is  a 

statement  of  the  fact  that  the  firequency  of  the  sinewave  must  be  less  than  or  equal  to  the 
Nyquist  frequency  determined  by  the  sampling  interval.  These  requirements  are  satisfied 
\i  At  = I IT  , N = 2^  and  k » j for  integer  values  of  j and  k ; for  these  choices  n . 

Generally,  the  best  accuracy  is  obtained  when  n = N 14  (j  = 2). 

As  an  example  of  this  technique,  the  sinewave 

xit,k)  = 1.25  sin  [2ji  (128.0)  r + (7.16) 


was  created  using  the  SINE  routine.  The  file  contained  ten  records  with  8192  points  per 
record.  The  sampling  interval  was  chosen  to  be  (using  n =N  14) 


At 


— — — = 0.001953125 
(4)(128) 


(7.17) 


so  that  the  record  contained  exactly  2048  cycles  of  the  sinewave.  The  spectrum  should  have 
a peak  at  a fi’equency  of  128  Hz  where  the  amplitude  is  given  by 


^^128)  = lim  E 

7-.® 


lim  X{f,T) 
128 


(1.25)(8192) 

(5 12) (2) 


(7.18) 


The  actual  value  reported  by  the  spectrum  routine  was  9.9656  at  a firequency  of  128.0082  Hz 
which  was  the  closest  discrete  frequency  value.  A plot  of  the  amplitude  spectrum  is  shown 
in  figure  7.8.  When  computing  this  spectrum,  no  mean  removal  or  detrending  was  used, 
and  no  filtering  was  performed.  Furthermore,  application  of  a tapering  function  would 
significantly  change  the  expected  amplitude  spectrum  and  therefore  no  tapering  was  used. 


The  required  inputs  to  the  SINE  routine  follow. 


sine 

Enter  the  amplitude  of  the  sinewave  : 1.25 
Enter  the  frequency  of  the  sinewave  : 128.0 
Enter  the  phase  (in  degrees) : 32.8 
Enter  DC  value  : 0.0 

Enter  delta  T secs.  (1.0/Sample  Rate)  : 0.001953125 
Enter  # of  points  per  record  (N)  : 8192 
Enter  # of  records  (ASIZE) : 10 
Enter  output  file  name  (4  chars) : SINE 


(gives  a simple  answer) 

(can  be  any  desired  number) 
(choose  any  initial  value) 
(no  DC  for  this  example) 
(according  to  eq.  7.17) 

(must  be  power  of  2) 

(OK  for  this  example) 


SINEWAVE  CREATION  UTILITY 


Creating  SINE.DAT  now. 

Writing  record  0010 

Use  voltage  transformation  factor  = 20.0  / 4096.0 
Program  terminated  successfully. 


60 


The  required  inputs  to  the  spectrum  routine  to  create  the  amplitude  spectrum  given  in 
figure  7.8  follow. 

spectrum  sine 

Transform  into  (V)oltages  or  into  (O)ther  quantity 
or  read  (R)eal  data  which  is  already  transformed  : v 

Enter  VOFST  1)  10/4096,  2)  20/4096,  3)  2/4096  or  4)  other  : 2 

Remove  (M)ean  value,  (D)etrend  or  (N)either  : n 

Filter  the  data  using  (L)owpass,  (H)ighpass 

(B)andpass,  Band(S)top  or  (N)one  : n 

How  many  channels  of  data  (1  or  2)  : 1 

(A)mplitude,  (P)ower  or  (Qross  spectrum  : a 

Taper  the  data  (Y/N)  : n 


Processing  SINE.DAT  now. 

Record  0010  : Ave  = .00000  Rms  = .88208 
File  Totals  : Ave  = .00000  Rms  = .88090 

No  tapering  applied. 

0010  records  will  be  used. 

Transforming  records  0009  and  0010. 

Writing  SINEA.PRN  now. 

starting  time  : 20:24:53.77 

ending  time  : 20:25:13.21 

elapsed  time  : 00:00:19.44 

Program  terminated  successfully. 

Note  that  the  maximum  number  of  points  per  record  is  controlled  by  the  value  of 
NMAX.  The  value  of  NMAX  is  currently  set  to  16384.  (If  the  user  should  ever  desire  to 
make  a change  to  the  value  of  NMAX,  the  change  must  be  made  both  in  the  SINE  routine 
and  also  in  the  SPECTRUM  routine.  NMAX  must  be  the  same  number  in  both  routines  at 
all  times.) 
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73  Experimental  Noise  Input 

The  input  data  file  for  this  example  was  derived  from  an  experiment  employing  an 
electrodynamic  shaker  table.  This  device  is  essentially  a large  solenoid  mated  to  an  electronic 
feedback  system  which  is  used  to  accurately  control  the  oscillatory  motion  of  the  core.  The 
displacement  amplitude  and  frequency  can  be  set  to  desired  levels  and  then  maintained  very 
accurately.  This  instrument  was  used  in  an  experiment  designed  to  test  the  frequency 
response  of  an  optical  system  which  was  built  to  measure  displacement  of  a surface  at  a 
point.  Prior  to  the  actual  frequency  response  test,  the  optical  system  was  used  to  measure 
the  displacement  of  the  core  of  the  solenoid  with  the  shaker  table  energized  but  undergoing 
no  motion.  The  computation  of  a power  spectrum  of  this  displacement  noise  can  then  be  a 
useful  diagnostic  tool  for  the  identification  of  noise  sources  in  subsequent  displacement 
spectrum  measurements. 

The  signal  was  digitized  using  a 12-bit  A/D  converter  with  an  input  range  of  -10  to 
10  V.  The  sampling  rate  was  2 kHz  which  resulted  in  a sampling  interval  of  At  = 0.0005 
seconds.  80  records  of  data  were  acquired  with  2048  points  per  record.  The  power  spectrum 
was  computed  for  a range  of  frequencies  up  to  the  Nyquist  frequency  which  was  1 kHz;  the 
resolution  was  given  by  A/=  UN  At  = 0.97656  Hz  . All  power  spectra  shown  in  this  section 
were  computed  by  first  detrending,  then  tapering  the  input  time  records  with  the  Welch  taper 
and  employing  50%  overlap.  The  resulting  power  spectrum  is  shown  in  figure  7.9;  loga- 
rithmic scales  are  used  for  greater  clarity.  The  figure  reveals  that  most  of  the  power  in  the 
signal  lies  below  100  Hz.  A 60  Hz  contribution  attributable  to  the  power  supply  may  be 
identified  as  well  as  other  resonances  and  their  harmonics. 
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This  input  file  allows  a useful  demonstration  of  the  filtering  capabilities  of  the  spectrum 
routine.  Shown  in  Figure  7.10  are  power  spectra  of  the  same  input  file  for  which  figure  7.9 
was  obtained;  however,  prior  to  the  spectral  calculation,  each  of  the  four  types  of  filters 
available  for  use  were  applied  to  the  data.  In  each  case  five  cascaded  filter  stages  were  used 
providing  a tenth  order  digital  Butterworth  filter.  As  can  be  seen,  the  rolloff  is  extremely 
rapid  with  power  in  the  stopband  3 to  4 decades  below  that  in  the  passband.  This  is  an 
effective  means  for  removing  undesirable  frequency  components  from  the  data  prior  to 
spectral  computation.  This  can  also  be  useful  for  quickly  estimating  the  amount  of  power 
contributed  by  specific  peaks  or  regions  of  the  spectrum.  The  rms  value,  calculated  as  the 
square  root  of  the  area  under  the  spectral  density  curve  (according  to  eq.  4.13),  is  given  in 
table  7.2  for  each  of  the  filter  cases  shown  in  figure  7.10  and  for  the  no-filter  case  plotted 
in  figure  7.9. 


No  filter 

Lowpass 

Highpass 

Bandpass 

Bandstop 

rms  values 

0.014803 

0.013319 

0.009529 

0.005894 

0.010116 

Table  7.2  RMS  values  reported  after  application  of  each  filter  type 

The  table  indicates  that  about  10%  of  the  contribution  to  the  root-mean-square  value  of  the 
signal  may  be  found  above  100  Hz,  about  a third  lies  below  30  Hz  and  just  under  half  may 
be  found  between  50  and  100  Hz.  However,  adding  the  rms  value  for  the  bandpass  case  to 
the  rms  value  for  the  bandstop  case  does  not  equal  the  rms  value  for  the  no  filter  case.  This 
discrepancy  is  due  to  the  fact  that  some  power  is  in  the  stopband  for  the  bandpass  and 
bandstop  cases,  and  this  contributes  a small  amount  to  the  rms  values  indicated.  Therefore, 
this  method  should  only  be  used  for  rough  estimates;  the  most  accurate  method  is  to  integrate 
the  spectral  density  function  over  fi-equency  bands  of  interest. 
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Figure  7.10  Effects  of  filtering  prior  to  spectrum  computation 
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7.4  Sampling  Considerations 

A full  discussion  of  sampling  procedures  will  not  be  reiterated  here;  a complete  dis- 
cussion of  these  techniques  can  be  found  in  the  references,  and  the  interested  reader  is  urged 
to  look  there.  Instead,  this  section  attempts  to  present  some  of  the  considerations  necessary 
when  making  choices  among  the  various  sampling  parameters,  and  the  manner  in  which 
these  choices  affect  the  spectral  density  results. 

The  most  important  parameter  is  probably  the  choice  of  sampling  interval  At . In 
general,  one  attempts  to  choose  this  quantity  to  be  of  the  order  of  the  smallest  time  scale  of 
interest  in  the  physical  phenomenon  under  investigation.  The  selection  of  At  determines 
the  largest  frequency  (bandwidth)  of  the  resulting  spectral  analysis.  This  critical  frequency 
is  known  as  the  Nyquist  frequency  and  is  given  by 


fc  = 


1 

2At 


(7.19) 


Another  quick  way  for  determining  the  Nyquist  frequency  is  to  note  that  fc^fsl'^  where 

^ = 1/At  is  the  sampling  rate  at  which  data  are  acquired.  If  energy  exists  in  the  signal 
above  the  Nyquist  frequency  determined  by  the  desired  sampling  rate,  then  this  energy  must 
be  removed  by  means  of  an  analog  lowpass  filter  prior  to  A/D  conversion.  This  is  necessary 
to  prevent  aliasing  of  the  high  frequency  components  in  the  data.  Aliasing  is  essentially 
confusion  among  the  high  and  low  frequency  components  in  the  data  that  occurs  whenever 
the  sampling  rate  is  set  too  low.  A very  good  discussion  of  aliasing  may  be  found  in  Bendat 
and  Piersol  (1986). 

Another  parameter  of  interest  is  the  number  of  data  points  per  record  per  channel  N . 
This  quantity  combined  with  the  choice  of  At  determines  the  record  length  T = N/St  and 
the  resolution  bandwidth  of  the  spectral  analysis. 


A/  = 


1 

NAt 


(7.20) 


A large  value  of  N for  a given  value  of  At  gives  a smaller  A/  and  can  give  a more  detailed 
spectral  density  estimate.  On  the  other  hand,  if  the  record  length  T is  very  long,  rea//zaZ?///ty 
problems  may  occur:  is  the  physical  phenomenon  stationary  for  these  long  time  periods? 
Is  there  a limitation  on  maximum  file  size?  Can  the  resulting  massive  amounts  of  data  be 
adequately  stored? 

Another  important  parameter  is  the  number  of  records  n^  of  length  T to  sample  for 
each  data  file.  These  records  should  be  contiguous  in  time  (in  order  to  use  overlapping)  and 
constitute  an  ensemble  of  statistically  stationary  records;  the  expected  value  operation 
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contained  in  the  definitions  of  the  spectral  density  estimates  is  approximated  when  an 
ensemble  average  is  calculated  for  this  set.  Recall,  that  it  is  (the  number  of  records 
before  overlapping)  which  is  used  in  the  determination  of  the  errors  associated  with  each 
of  the  spectral  density  estimates  and  not  N ot  ^ . Increasing  rapidly  decreases  the 
normalized  random  error  e up  to  about  = 100  or  so;  for  larger  values  of  the  error 
decreases  much  more  slowly.  There  may  be  certain  situations  in  which,  for  a given  overall 
file  size,  it  is  better  to  increase  and  decrease  N for  a given  A/  in  order  to  increase 
accuracy  at  the  expense  of  resolution  bandwidth  A/ . 
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8 Utility  Routines 

Four  utility  programs  are  included  on  the  diskette.  The  SINE  and  COLOR  routines 
have  been  discussed  in  the  previous  section  concerning  examples.  MAKDAT  is  a routine 
designed  to  read  a data  file  in  ASCII  format  and  then  to  write  the  necessary  file  header  and 
data  to  an  unformatted  file.  MAKCOEF  is  similar;  it  reads  a user-prepared  ASCII  format 
coefficient  data  file  and  then  writes  the  file  header  and  coefficients  to  an  unformatted  file. 
It  is  easiest,  of  course,  to  add  the  necessary  code  to  create  the  unformatted  files  directly  to 
the  routine  that  is  used  to  create  (or  collect)  the  data  in  the  first  place;  however,  if  the  user 
is  unfamiliar  with  unformatted  FORTRAN  data  files  or  the  FORTRAN  language  in  general, 
then  these  routines  attempt  to  solve  the  problem. 

8.1  MAKDAT  routine 

In  order  to  use  MAKDAT,  the  data  to  be  analyzed  by  the  spectrum  routine  must  be 
generated,  or  collected  by  an  A/D  converter,  and  then  stored  to  a file  in  ASCII  format.  In 
this  format  the  file  may  be  viewed  by  any  ASCII  editor,  can  be  displayed  on  the  screen  and 
can  be  printed  on  a printer.  This  is  not  an  optimum  way  to  store  data;  these  files  tend  to  be 
extremely  large.  In  anticipation  of  this  fact  other  types  of  data  files  were  explored  for 
possible  use,  and  unformatted  FORTRAN  files  appeared  to  be  the  simplest  solution. 
Unformatted  files  are  typically  a factor  of  five  smaller  in  size  than  the  equivalent  ASCII 
file.  An  example  will  be  given  to  illustrate  the  means  for  transferring  the  data  to  unformatted 
files  which  can  then  be  used  by  SPECTRUM. 

Suppose  that  100  records  of  one  channel  of  data  saved  as  2-byte  integers  and  sampled 
at  10000  samples  per  second  with  2048  data  points  per  record  have  been  acquired,  lliese 
data  should  be  placed  in  an  ASCII  file  with  the  file  extension  [.ASC]  and  with  a file  name 
conforming  to  the  various  naming  conventions  discussed  earlier.  The  file  must  contain 
numbers  only;  the  numbers  may  be  written  one  per  line  or  many  per  line  (but  this  does  not 
shorten  the  length  of  the  file).  Blank  lines  may  be  added  if  desired.  The  first  number  in  the 
file  should  be  the  first  data  point  of  the  first  record;  the  2049th  number  in  the  file  should  be 
the  first  data  point  of  the  second  record,  and  the  last  number  in  the  file  should  be  the  2048th 
data  point  of  the  100th  record.  A session  with  MAKDAT  to  accomplish  this  transformation 
would  proceed  as  follows. 

makdat 

(I)nteger  (2-byte)  or  (F)loating-point  (4-byte)  data  : i 

Enter  # of  channels  (1  or  2) : 1 

One  channel  : Total  # points  per  record  <=  16384. 

Two  channels  : Total  # points  per  record  per  channel  <=  8192. 

Enter  # of  points  per  record  (power  of  two)  : 2048 

One  channel  : Must  be  EVEN  # of  records. 
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Two  channels  : May  be  EVEN  or  ODD  # of  records. 

Enter  # of  records  in  the  data  file  : 100 

One  chan  : Delta  t is  spacing  between  data  points. 

Two  chans  : Delta  t is  spacing  between  data  pts  - SAME  channel 
Delta  t divided  by  2 is  spacing  between  data  pts 
- different  channels. 

Enter  sampling  interval  delta  t (secs)  : 0.0001 
Enter  ASCII  input  file  name  (4  chars)  : test 

DATA  FILE  CREATION  UTILITY 
Creating  TEST.DAT  now. 

# channels  = 1 

record  size  = 4096  bytes 

# of  records  =100 

delta  t = 100  microseconds 

Record  0100 

Program  terminated  successfully. 

The  MAKDAT  routine  reads  the  user-prepared  TEST.ASC  file  and  calculates  the  values 
needed  for  the  file  header.  The  file  header  would  consist  of 

1 4096  100  100 

The  routine  then  transfers  the  data  to  the  unformatted  file  with  the  name  TEST.DAT.  Finally, 
to  process  this  file  with  the  spectrum  routine  the  user  would  enter:  spectrum  test. 

For  two  channels  of  data,  it  is  necessary  that  the  numbers  be  stored  in  each  record  of 
the  ASCn  file  in  the  following  manner: 

one  channel:  XqX^X2  ... 

two  channels:  ^oyo^iTi^>’2  ••• 

using  the  notation  defined  earlier.  Therefore,  the  first  number  in  the  file  will  be  the  first 
data  point  of  the  first  record  for  channel  1;  the  second  number  will  be  the  first  data  point  of 
the  first  record  for  channel  2;  the  third  number  in  the  file  will  be  the  second  number  for  the 
first  record  of  channel  1,  and  so  on.  Also,  care  should  be  taken  to  enter  correct  responses 
to  the  prompts  when  two  channels  of  data  are  to  be  transferred.  In  response  to  the  number 
of  points  per  record  per  channel  prompt,  the  user  would  enter  1024  points  if  the  data  file 
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consisted  of  two  channels  of  data  with  2048  points  per  record  and  therefore  1024  points  for 
each  channel  in  each  record.  Note  that  the  maximum  number  of  points  per  record  is  controlled 
by  the  value  of  NMAX;  the  maximum  number  of  points  per  record  per  channel  is  then 
NMAX  / 2.  The  value  of  NMAX  is  currently  set  to  16384.  (If  the  user  should  ever  desire 
to  make  a change  to  the  value  of  NMAX,  the  change  must  be  made  both  in  the  MAKDAT 
routine  and  also  in  the  SPECTRUM  routine.  NMAX  must  be  the  same  number  in  both 
routines  at  all  times.)  Similarly,  when  responding  to  the  prompt  for  the  sampling  interval, 
the  user  would  enter  0.0002  seconds  for  the  sampling  interval  between  data  points  of  the 
same  channel  for  the  case  in  which  the  data  were  sampled  at  10000  samples  per  second. 
Although,  0.0001  seconds  is  the  sampling  interval  between  adjacent  data  points  in  the  record, 
the  sampling  interval  for  data  points  of  the  same  channel  is  a number  twice  as  large. 

8.2  MAKCOEF  routine 

The  primary  reason  for  using  unformatted  coefficient  data  files  with  the  spectrum 
routine  is  to  allow  the  user  a means  for  preparing  all  of  the  coefficient  data  for  the  trans- 
formation of  voltages  into  other  dimensional  quantities  for  several  files  which  will  then  be 
processed  by  SPECTRUM  in  a batch.  The  MAKCOEF  utility  allows  the  user  to  prepare 
this  coefficient  data  in  an  ASCII  file;  the  utility  will  write  the  file  header  and  then  transfer 
the  data  to  an  unformatted  coefficient  data  file  which  can  be  read  by  the  spectrum  routine. 

Suppose  that  it  was  desired  to  process  the  data  files  D061.DAT  thru  D086.DAT  in  a 
batch  by  the  spectrum  routine.  Each  of  these  data  files  consisted  of  two  channels  of  2-byte 
integers  taken  fi’om  an  experiment  in  fluid  mechanics,  say.  Furthermore,  suppose  that  the 
quantities  sampled  were  two  components  of  a fluid  velocity,  which  were  transduced  into 
voltages  by  some  means  and  then  sampled  by  an  A/D  converter.  The  user  knows  the  cal- 
ibration data  for  the  transducer,  and  the  required  transformation  fi’om  voltages  back  into 
velocities  can  be  accomplished  for  each  velocity  component  by  using  a polynomial  of  up 
to  fourth  degree  in  the  form: 


+ b^vl  + b^vl  + by^ 


(8.1) 


where  v„  is  the  nth  data  point  expressed  as  a voltage,  bQ,b^,...,b^  are  the  coefficients  of 

the  velocity  transformation  and  u„  is  the  nth  data  point  which  has  been  converted  to  a 
velocity.  Now,  suppose  further  that  the  data  in  each  of  the  files  were  taken  over  a period 
of  time  in  which  the  calibration  coefficients  changed;  thus,  a different  set  of  coefficients  are 
required  for  each  of  the  26  files  to  be  processed.  Following  the  instructions  in  the  section 
on  coefficient  data  files,  this  scenario  would  require  the  creation  of  two  unformatted  coef- 
ficient data  files,  DCON.DAT  and  DCON2.DAT,  each  containing  the  26  sets  of  five 
coefficients  bQ,bi,...,b^  for  the  velocity  transformation  for  channel  1 and  channel  2, 
respectively. 
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The  user  should  create  each  ASCII  data  file  with  a four  character  name  followed  by 
the  extension  [.ASC].  The  names  should  conform  to  the  various  naming  conventions 
specified  earlier.  The  ASCII  file  must  contain  numbers  only,  although  blank  lines  may  be 
inserted  into  the  file  if  desired.  The  numbers  will  be  transferred  to  4-byte  floating  point 
quantities  in  the  unformatted  data  file  thereby  limiting  the  coefficients  to  single  precision 
(6  or  7 decimal  places).  The  coefficients  may  be  stored  one  per  line  or  many  per  line  in  the 
ASCn  file.  The  first  number  in  the  file  will  be  bo  for  channel  1 of  data  file  D061.DAT; 
the  fifth  number  in  the  file  should  be  for  channel  1 of  data  file  D061.DAT;  the  sixth 
number  in  the  file  will  be  for  channel  1 of  data  file  D062.DAT,  and  so  on.  The  130th 
and  last  number  in  the  file  should  be  b^  for  channel  1 of  data  file  D086.DAT.  The  user 
should  then  create  a second  ASCII  file  in  a similar  manner  specifying  all  of  the  coefficients 
for  channel  2 of  the  26  files  to  be  processed.  The  utility  MAKCOEF  should  then  be  executed 
twice,  first  transferring  the  contents  of  the  first  ASCII  file  (for  channel  1)  to  DCON.DAT 
and  then  for  the  second  ASCII  file  (for  channel  2)  to  DCON2.DAT.  A session  using 
MAKCOEF  would  proceed  as  follows. 

makcoef 

Enter  first  letter  of  data  file  to 

which  these  coefficients  will  be  associated  : d 

Are  these  coefficients  for  channel  (1  or  2)  : 1 

Enter  # of  sets  of  coefficients  : 26 

Enter  starting  file  number  : 61 

Enter  ASCII  input  file  name  (4  chars)  : chnl 

COEFFICIENT  FILE  CREATION  UTILITY 

Creating  DCON.DAT  now. 

# sets  of  coeffs  = 26 

# coeffs  in  each  set  = 5 

# of  starting  file  = 61 

Program  terminated  successfully. 

The  MAKCOEF  routine  will  read  the  ASCII  file  CHNl.ASC  and  will  write  the  following 
file  header  to  the  unformatted  file  DCON.DAT: 


26  61 
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The  coefficients  for  channel  1 are  then  transferred  to  the  unformatted  file.  Note  that  5 
coefficients  per  set  are  required;  therefore,  this  number  is  not  written  to  the  file  header.  If, 
at  a later  time,  the  user  wishes  to  reprocess  the  files  D074.DAT  through  D080.DAT  (perhaps 
specifying  different  options  in  the  spectrum  routine),  it  is  not  necessary  to  re-create  the  two 
coefficient  data  files;  as  long  as  the  data  files  to  be  reprocessed  are  chosen  from  among  the 
26  original  members  of  the  set,  the  spectrum  routine  will  use  only  the  coefficients  that  it 
needs.  Similarly,  if  the  user  wishes  to  process  the  data  files  in  an  arbitrarily  numbered  order 
(not  consecutively),  the  coefficient  data  files  do  not  need  to  be  changed  as  long  as  the  data 
files  come  from  the  26  original  members  of  the  set.  Thus,  D061.DAT,  D064.DAT, 
D069.DAT  and  D086.DAT  could  be  processed  by  SPECTRUM  without  changing  the 
coefficient  data  files.  The  converse  is  also  true,  however.  If  the  user  desires  at  the  outset 
to  process  only  the  four  arbitrarily  ordered  files  above,  two  coefficient  data  files  must  be 
created  which  contain  not  4,  but  26  sets  of  5 coefficients  for  the  data  files  D061.DAT  through 
D086.DAT.  The  spectrum  routine  requires  that  sets  of  coefficients  in  the  coefficient  data 
file  be  listed  for  consecutively  numbered  files  even  if  not  all  of  the  sets  of  coefficients  are 
used.  With  this  in  mind,  it  is  perhaps  best  to  create  the  coefficient  data  files  for  an  entire 
set  of  input  data  files;  then  members  of  the  set  can  be  processed  by  the  spectrum  routine 
either  consecutively  or  in  an  arbitrarily  numbered  order  at  any  subsequent  time.  Of  course, 
the  user  could  transform  the  voltages  back  into  dimensional  quantities  using  his/her  own 
code  and  store  the  resulting  floating-point  numbers  into  ASCII  data  files  (then  use  MAK- 
DAT)  or  directly  into  unformatted  data  files  to  be  processed  by  SPECTRUM;  then,  all  of 
this  mess  concerning  coefficient  data  files  could  be  avoided  altogether. 

Note  that  the  maximum  number  of  points  per  record  is  controlled  by  the  value  of 
NMAX.  The  value  of  NMAX  is  currently  set  to  16384.  The  maximum  number  of  files  that 
may  be  processed  by  SPECTRUM  in  a batch  is  controlled  by  the  setting  of  IFMAX;  this  is 
currently  set  to  100.  (If  the  user  should  ever  desire  to  make  a change  to  the  value  of  NMAX 
or  IFMAX,  the  change  must  be  made  both  in  the  MAKCOEF  routine  and  also  in  the 
SPECTRUM  routine.  NMAX  and  IFMAX  must  be  the  same  numbers  in  both  routines  at 
all  times.) 
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